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ABSTRACT 


The Black-Scholes model is commonly used to track the price of European options 
with respect to maturity in many financial markets. This model degenerates into 
a partial differential equation that relates the European-style option price to the 
underlying price and time of expiry. Black-Scholes assumes that underlying prices 
satisfy a geometric Brownian motion. 

After the U.S. stock market crash of 1987, this assumption becomes inaccurate as 
it fails to represent the behavior of S&P 500 European vanilla option prices. Specifi- 
cally, under the measure of moneyness, the volatility smirk does not flatten out and 
the resulting conditional return distribution does not converge to normality. Recent 
academic literature have proposed readjusted financial models to account for the 
shortcomings of Black-Scholes, none which successfully have combined infinite return 
moments and finite price moments. 

To reduce the effects of these consequences and to incorporate the additional 
moment conditions, we assume that the underlying satisfy a Levy a—stable motion. 
Under this assumption, we will derive the Finite Moment Log Stable (FMLS) model 
and its respective fractional partial differential equation counterpart. Then, we will 
solve the Black-Scholes equation under FMLS by using the standard finite difference 
method and a finite volume scheme that significantly reduces the computational and 
storage cost in comparison. Lastly, we will perform a numerical simulation of our 
methods by using recent financial data in the S&P 500 market acquired within a 


one-year time frame to compare the performance of these methods. 
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INTRODUCTION 


Constructed by Fischer Black and Myron Scholes, the Black-Scholes model is a fi- 
nancial instrument designed to analyze the behavior of European option prices under 
changes of maturity dates and the stochastic behavior of asset prices. Contrary to 
most financial markets, this model assumes that asset prices are distributed normally 
and thus, satisfy a geometric Brownian motion. Realistically, the volatility and in- 
terest rate are varying and randomly behaved. To simplify the construction of the 
Black-Scholes model, we assume that these parameters are constant and known. Our 
choice of these parameters will be highly dependent on recent trends of the current 
market. [2] Assuming the conditions previously mentioned, there exists a partial dif- 
ferential equation that captures the behavior of European option prices satisfying the 
Black-Scholes model. In most financial and mathematics literature, this equation is 
commonly known as the Black-Scholes equation. The derivation of Black-Scholes can 
be reproduced with stochastic analysis using applications of finance. One can also 
derive the closed form solution of the Black-Scholes model by solving the traditional 
diffusion model under specific changes of variables and Fourier transforms. 

However, the Black-Scholes model fails to capture the behavior of option prices 
under rare and extreme circumstances as these assumptions are not applicable real- 
istically. This is shown immediately after the U.S. stock market crash of 1987, where 
numerous academic resources have identified a consistent pattern found in the S&P 
500 options market. At a given maturity date, the Black-Scholes model implies that, 
with respect to the strike price, the volatilities for out-of-the-money puts are much 


higher than out-of-the-money calls. In finance, this phenomenon is known as volatil- 


ity smirk. Under the measure of moneyness, the implied volatility smirk does not 
flatten out as maturity increases. Namely, one can observe that the downward slope 
of the volatility smirk corresponds to the asymmetry of the return distribution and 
the positive curvature of the smirk correlates with the existence of leptokurtosis. [3] 
Therefore, the conditional return distribution does not converge to normality, the 
central limit theorem does not hold, and the assumption that asset prices satisfy a 
geometric Brownian motion cannot be made in these rare situations. This calls for 
a different model that accurately tracks this recent behavior of the S&P 500 options 
market. 

In doing so, we require the construction of a variant of the Black-Scholes model 
to generalize the behavior of option pricing in various financial markets. We need to 
model returns by considering an a-stable motion with maximum negative skewness 
where aq is the tail index satisfying 0 < a < 2. The result of this model will combine 
infinite return moments and finite price moments. These adjustments are necessary 
to accurately reflect the observed behavior of S&P 500 option prices documented after 
the market crash. When a = 2, the a-stable motion becomes the standard Brownian 
motion, and the resulting model constructed from the a-stable motion degenerates 
into the standard Black-Scholes model. Setting a < 2 yields many interesting prop- 
erties that lie in the foundation of the study of fractional calculus. Similarly to the 
Black-Scholes equation, there exists a fractional counterpart that reflects the generic 
behavior of European vanilla option pricing. We begin Chapter 1 by deriving the clas- 
sical Black-Scholes equation using methods in stochastic calculus and its applications 
in finance. Then, we will gather an intuitive understanding of fractional derivatives 
to study the traditional diffusion model and its fractional counterpart. We will con- 
struct a generalization of the Black-Scholes model that requires the theory of stable 
random variables to accurately track the behavior of European option pricing. Un- 


der this derived model, we will finally derive and analyze the Black-Scholes equation 


under this model. 

After this derivation, observing the behavior of the price of a European vanilla 
option under this model and determining the respective analytic solution becomes 
a natural question. In doing this, we study the closed form of the Levy density 
function, which can be written in terms of the Fox H-function. When a = 2, the Levy 
density function becomes the probability density function of the standard Gaussian 
distribution. This function plays a crucial role in determining the analytic solution of 
the price of European vanilla options. The properties that arise from the analytic form 
of the Fractional Black-Scholes equation and its consistencies with the Black-Scholes 
formula become other questions of interest. In Chapter 2, we derive this analytic 
form using the Fox H-function and propose associated theorems that emphasize the 
properties of the closed form corresponding to the Black-Scholes formula. We prove 
these theorems to provide mathematical justification for the observations documented 
in many options markets. We will also analyze the asymptotic behavior of the closed 
form solution with respect to the log price and show that this solution satisfies the 
put-call parity. 

Many numerical methods have been proposed to solve partial differential equa- 
tions with fractional order. Among many, the finite difference method is known to 
be an effective method that uses numerical values evaluated at specific nodes to ap- 
proximate intermediate values. This is achieved by partitioning bounded intervals 
into finite increments and using the finite difference approximation of the derivatives 
to construct a recursion that defines the behavior of the function. To use the finite 
difference method, we are required to truncate the unbounded interval to the finite 
interval |—y, y] for some constant y. The constant y depends on the maximum price 
of the asset and can be readjusted with respect to the necessary conditions of the 
market. In Chapter 3, we recall the definition of a fractional derivative in terms 


of Grunwald weights from Chapter 1. We will show that a shift on the weights is 


required to yield a consistent and unconditionally stable solution of the finite differ- 
ence method. Under the assumption that the log price is contained in some bounded 
interval, this will result in a recursion that relates numerical values of consecutive 
nodes. This recursion can be reverted into a matrix equation which can be solved 
using standard Gaussian elimination. 

Determining more efficient methods has been an area of interest in the study of 
numerical analysis. We propose a few approaches in using efficient finite difference 
methods which reduces the computational cost from O(m*) to O(m log? m) and the 
storage cost from O(m?) to O(m) where m represents the number of unknowns. [13] 
Contrary to the standard method of Gaussian elimination, this reduction of cost can 
be achieved by proposing banded coefficient matrices that approximate the default 
coefficient matrix with reduced storage cost. In the first approach, we require the 
boundary conditions to take the value of zero at the endpoints. We circumvent this 
issue by using the Taylor series expansion of the closed form to determine two points 
which will take the same values. In the second method, we require the partial dif- 
ferential equation to be rewritten conservatively. Both methods will generate matrix 
equations that can be solved using banded coefficient matrices. In Chapter 4, we will 
first manipulate the current problem so that the necessary conditions for both meth- 
ods hold. After doing this, we will apply these methods for the fractional counterpart 
of the Black-Scholes equation to determine the behavior of European vanilla option 
pricing under the derived Finite Moment Log Stable process defined in Chapter 2. 
Lastly, we will implement a Fast Conjugate Gradient Squared Method to accelerate 
the performance of the finite volume scheme by reducing the storage cost. 

In the last chapter, we will perform a numerical simulation of the finite difference 
methods discussed throughout this thesis by using data observed in the most recent 
one year time period of the S&P 500 options market. We will simulate the prices 


in a graph to observe the effects of European put options as a@ varies. We will also 


discuss the accuracy of the finite volume scheme and emphasize the low CPU time 
required to implement the method under various approaches. We will finally make 
observations on the data acquired in our simulation to support the theory developed 


in the previous chapters. 


CHAPTER 1 


MODEL PROBLEM 


In this chapter, we will derive the Black-Scholes equation using financial applications 
and stochastic calculus. This will require many assumptions to simplify the derivation 
of the model. Under these assumptions, we will discuss why the Black-Scholes model 
fails to capture the recent behavior of the S&P options market. We will proceed to 
acquire an understanding of fractional derivatives by analyzing the fractional coun- 
terparts of the traditional diffusion model and the central limit theorem. We turn to 
a generalized class of stochastic processes that mimics many similar characteristics 
of the conditional return distribution in the S&P options market. Specifically, we 
will define and study properties of a-stable motion and derive a readjusted model 
highlighting the inconsistencies of the traditional Black-Scholes model. Then we will 
conclude this chapter by deriving the Black-Scholes equation under the FMLS model 


using our understanding of fractional calculus and a-stable motion. 


1.1 PRELIMINARIES 


In 1826-1827, Robert Brown conducted a study on the movement of particles from 
pollen grains suspended in water. He noted that these particles move in an irregular 
manner and that the motions of two distinct particles appear to be independent. 
These observations are similarly observed in fluctuations of stock prices by French 
mathematician Louis Bachelier in the early 1900s. 

In 1905, Albert Einstein provided a mathematical explanation of the phenomena. 


Consider a tube of clear water and inject a unit amount of ink at time t = 0 and 


at location x = 0. Define u(x,t) to be the density of ink particles at location x and 
time t. Let f(y,7) be the probability density function of an ink particle moving from 
distance x to x + y in relatively small time rt. Then, by Taylor Series expansion of 


u(x — y,t), we have: 


(6S cebu") Flusridy 
= fe (SP) paride + (70 


ee Noe ax” 
lee) 0 fag 
=f (u(x,t) — Mey + 5 soy f(y,7)dy + O(r) 
leo) 0 leo) 
= ule.t) [ fy.rdy 29 Py py, a)dy 


Since f is a probability density function, for any small 7, we have / faajdy= 


1. Given f(—y,7) = f(y,7), it follows that / yf (y,T)dy = 0. We assume that the 


variance of f is linear in T; that is, for some constant o? > 0, we have, 


i y f(y, 7)dy = 07. 
Continuing from our calculations, we have, 


u(o,t+7) =ule.t) [ flysridy— M29 I y py, ray 


LOPE ey) of 05 
a 5 oe fo f(y, 7)dy + O(7) 


LO ua to 
5 Bx2 or + O(r) 


= u(x,t) 4 


Equivalently, 


= 2 
u(x,t +7) — u(z, t) - 10 u(x,t) + O(r). 
T 2 O77 


Letting tT > 0 yields, 


u(x,t) 07 O?u(z,t) 


Ot 2 Ox? 


The above partial differential equation is commonly known as the traditional 


diffusion model. This model has the following closed-form solution: 


1 x? 
exp | ——— ]. 
2Qra7t ; 207t 


This diffusion model and its solution can be mathematically justified using random 


wnt). = 


walks and the Laplace-DeMoivre Theorem. A rigorous proof of the derivation and 
the solution of this model can be found in Appendix A. 

Later in this thesis, we will provide a second proof of the diffusion model using 
Fourier transforms. We can extend this alternate proof to create a diffusion equation 
involving fractional derivatives to gather intuition of super-diffusion. Specifically, we 
will generalize our understanding of diffusion to develop an understanding of fractional 
derivatives and their properties. This understanding is crucial to the derivation of 


the Black-Scholes equation under the FMLS model. 


1.2 DERIVATION AND CONSEQUENCES OF BLACK-SCHOLES MODEL 


In this section, we will derive and discuss the consequences of the Black-Scholes 
model. More information can be referred to Roberts. [11] Let C(.S,t) be the price of 


a European call option for asset price S' and time t. We assume that C(S,t) varies 

C OC OC 
an 

Ot? OS’ OS? 


smoothly with respect to t and S' so that its partial derivatives 


are well defined and smoothly varying. 

We wish to construct a risk-free portfolio of one call option and y units of assets. 
Let II be the value of this portfolio. Since I] = —C(S,t) + yS is a function of 
stochastic asset value S’, II is an Ito process. To observe the behavior of II over a 
small incremental change of time dt, we must consider the stochastic differential of 


II, namely dll = —dC'+ gdS. 


The key assumption of the Black-Scholes model is that asset price satisfies a 


geometric Brownian motion. Therefore, we have 
ds 
aes = adt + BdW; 


where W;, represents Brownian motion with respect to time t, and a and £ repre- 


sents the stock drift and stock volatility respectively. Equivalently, 
dS = aSdt + BSdW;. 


We will use Ito’s Chain Rule (Theorem B.1) to derive the Black-Scholes equation. 
The proof of this result can be found in Appendix B. Applying Ito’s Chain Rule with 
Y =-C,x2=S, w=aS, and o = BS yields the following: 


OC oC 1 OC OC 
Tl — _ _ _ 2a2 _ po 
d (( ae — aS — SBS a) dt asscaw) + y(aSdt + BSdW,) 
A 0e. - Ae tins Se, ac 


Under the assumption that the portfolio is risk-free, we require the volatility of 
dll to be zero. This implies that to guarantee this assumption, we need y = Th In 


doing so, we obtain the following: 


OC: . Ths sear 
ant = ( — 509s ) a 


Ot 3? 
Furthermore, the risk-free assumption implies that the returns of the portfolio 
equal to the returns of the investments made in bonds. With the value of the portfolio 
being —C + eae and interest rate r, we have 


as 
ac 
T=r(-C+——S) dt. 
d r( c+ 555) at 


Equating the coefficients and rearranging terms yields the renowned Black-Scholes 


equation for the value C(S,t) of the call option, 


With the boundary conditions that forO <t<Tand0<S<w,C(S,t)h~S 
as S — oo, and C(S,T) = max{S — K,0} for strike price K, the solution of the 


Black-Scholes equation, known as the Black-Scholes Formula, is provided below: 


C(S,t) = SN(d,) — Ke? N(da) 


jheré, 
NG = i oe (-53°) ds 
And, 
4, — loalS/K) + (r+ PVT -2) 
BVT ot 
iy BODEN og _ gy 


Remark: Note that: 


B-B_1 
2 2 


(dy dz)(d, — dz) => SBVT 8 (2d — BVT -t) = log(S/K) +r(T —t). 


Therefore, 


Sart) _ N'(d2) 


Ko -sN“(di) 


<=> SN'(d}) = Ke"? N' (dp). 


The derivation of the Black-Scholes formula requires the correct substitutions for 
many variables to transform the Black-Scholes model into the traditional diffusion 
model. More details on the derivation of the Black-Scholes formula can be found in 
Appendix C. 

The Black-Scholes equation emphasizes the relationship between many observable 
parameters and the stock volatility (6; specifically by using the Black-Scholes Model, 
6 can be determined by the observable call option value C(S,t), the asset price S, 


the time to maturity t, and the risk-free interest r. 


10 


Note that by fixing these observable parameters in the Black-Scholes formula, we 


have: 
OC Od, = Odz 
eee eS, N’ pean ee K r(T t) NV! ae 
Od, 


= VT -t (Ke) N'(d2)) + 36 (SIN'(d1) — Ke") N'(d)) 


By the above calculations, we can deduce that the value of the call option is a 
monotonic increasing function of its implied volatility. This enables us to compute the 
volatility for various call options with respect to different asset prices and maturities. 
Due to monotonicity, the implied volatility can be interpreted as a re-scaling of option 
prices necessary to analyze its behavior with respect to the observable parameters. It 
is known that in most financial markets, the implied volatility flattens out as maturity 
increases, and thus, the conditional return distribution should converge to normality. 
In finance, this behavior is renowned as the volatility smirk. 

However, trends following the 1987 U.S. stock market crash suggest otherwise. 
Contrary to the implications of the Black-Scholes model, it is well documented that 
according to the S&P 500 index options market, the implied volatility does not flatten 
out as maturity increases. This is made apparent when implied volatility is observed 
with respect to maturity and financial measure moneyness d := sas), for some 
constant C’. Moneyness approximately determines the number of standard deviations 
that the strike is apart from the forward price. It can be observed that the downward 
slope of the observed smirk corresponds to asymmetry and the positive curvature of 
this smirk corresponds to the fat tails in the conditional return distribution. Due to 
these observations, the Black-Scholes model fails to capture the true behavior of the 


recent trends of option pricing. 


11 


This motivates us to study the effects of the return distribution when one changes 
the stochastic behavior of asset prices. Recall that for the Black-Scholes model, we 


assumed the following: 


e The value of the call option varies smoothly. 


e The asset price satisfies a geometric Brownian motion and returns are lognor- 


mally distributed. 


e The portfolio is risk-free and the model assumes constant volatility. 


These assumptions are consequences of the validity of the central limit theorem. 
We wish to readjust our model to compensate for the shortcomings of Black-Scholes. 
In most financial literature, it has been proposed to assume that asset prices satisfy 
a more generic class of stochastic processes called Levy a—stable motion. More 
discussion on this class of stochastic processes will be presented later in this chapter. 
One can construct a variant of the Black-Scholes equation that accurately reflects the 
behavior of implied volatility over maturity in many financial markets. 

In the next section, we will reprove the diffusion model using Fourier transforms 
and highlight the correlation of this model to the central limit theorem. We will dis- 
cuss the proof conceptually by using the theory of random walks and extend this proof 
to derive the fractional diffusion model. Then, we will briefly discuss the correlation 


of the fractional diffusion model to the conceptual ideas of super-diffusion. 


1.3. TRADITIONAL AND FRACTIONAL DIFFUSION MODELS 


Recall that the traditional diffusion model is governed by the partial differential 
equation: 


Ou(z,t) 07 O?u(z,t) 


Ot 2. cor 


and has the closed-form solution: 
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1 ae 
ut) = re 353 | 


We will reconstruct a proof of the diffusion model using Fourier transforms. This 


will require a definition and an example that will be useful later in this thesis. 


Definition 1.1. Let Y be a random variable with probability density function f(x), 
differentiable and bounded on the interval [a,b] with f™(a) = f(b) for all n. We 


define the Fourier transform of f(x) to be the following: 


fk) = FUF@)} = Ble” = [ Y emikt (nd. 


Example 1.1. We will attempt to compute the Fourier transform of f)(x). We 


A 


claim that f(a) = (ik)"f(x) and proceed to prove the claim with induction. By 


definition and applying integration by parts on f’(x), we have, 
fin) = [ef @da 
= [eo f(alt + ik fe fla)ar 
= Ahy (ae 
Now assume f(x) = (ik)™f (x). Then, similarly, 


a 5 , 
flr) (k) = / e the £(M4+)) (4) der 
‘ : Bs 
= [ef (ay + ik [em F(a hd 
= (ik) f(z). 
Therefore, we can deduce that f(x) = (ik)" f(x) for all n. 


Let {Y,,} denote a sequence of normal independent and identically distributed 
random variables that represent the jumps of a randomly selected particle. We define 
a random walk as the random variable, S, = Y, + ---+ Y,, which represents 


the position of the particle after n steps. We can reinterpret this random walk by 


13 


considering the position of a particle at time t > 0 and some scaling constant c (i.e. 


let n = ct). Furthermore, suppose that E[Y,] = 0 and E[Y,?] = o? > 0 and let f(z) 


n 


be the probability density function of Y,,. By considering the Taylor series expansion, 


we can deduce the following: 


f(k) = fe™ Fade 
= y (1 —ikax + Seay +: .) f(x)dax 


2! 
= [foe — ik [xf (v)de = < [o?f(a)de + o(k”) 
aps si +o0(k") 


where the integrals are evaluated over the defined range of x. Note that by 


linearity of expectation, we have, 


Ele Sn] = Ble eater + Fn) = E[e~*** | oo Ele Yn] = Ble" = f(k)”. 


Let u(x,t) be the probability density function of a random particle at time t > 0 
and distance x. Since u is a probability density function, we require the scaling 


-1/2 to normalize S,. Therefore, 


oe 2f.2 n 
Ele Sn) = (1 ae + xe) ; 


constant c 


By using (1 Sie oe o(n-1))" — e’, letting c > oo and n > ow yields, 


exp (— 3107?) = Ele] = G(k,t) 


where Z is a normalized random variable. It is plain that a(k,t) solves the fol- 
lowing differential equation: 


2 


dii(k, t) ae 
2 


(ik)2a(k, t). 


2 
= -FPalk,t) = 
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Inverting the differential equation yields the traditional diffusion model. Similarly, 
inverting the solution of the differential equation yields the closed-form solution of 
the diffusion model. More details on the inversion are supplemented in Appendix D. 


Consequentially, note that by Theorem D.1, since 


og 1 ; ; i 1 
Ele tke Sn) — exp (- tok?) = Ele 7] = eke Soh exp (—5t07#*), 
TO 


we have, 


Vi doseasy, 
(ig. = IE 


Jn 
when Z is Brownian motion with mean zero and variance o*t. The convergence 
of the distribution of random variables is commonly renowned as the central limit 
theorem. 
We will now extend this proof to derive the fractional diffusion model. Let X be 
a Pareto random variable. Then, for all n, PLY < 2] = 1—Ca~ for constant C > 0, 
x > C/* and 1 <a < 2. Taking the derivative of both sides yields the probability 


density function denoted as: 


Cage? x > Cre 
f(x) = 
0 g< Ce 


Let 0 < p< a. Then the pth moment is given in the following computation: 


E[X?] = [efode 
= Ca oe Polder 
bas I. a 


Pp—a@ Cl/a a—p 


CP/o 


= ca| 


Given that 1 < a < 2, the first moment exists and the second moment is unde- 
fined. This implies that the mean of X exists and the variance of X doesn’t exist, 
and thus, the central limit theorem does not hold. We will have to take this into 
consideration as we construct the fractional diffusion model. 


We state the following proposition: 


15 


Proposition 1.1. Let X be a Pareto random variable with probability density function 
defined as above for some 1<a< 2. Then, as k + 0, the Fourier transform of f(x) 


1s 


e —1_qi/a ry 7 
fk) =1- C9 ik +0 


Proof. Note that, 


Ele **] = : et Ogg % dx 
cl a 


= a [1 —ikax + (ex —1+ ike) Cax~* dx 


=1-cVve = 


ik + i. (er —1+ ike) Cax~*"!dx 
a-—l 0 


Cil/e : 
— | (e7i# —1+ ike) Caz °"dzx. 

0 

Let 
J(a):= or i (or =f ike) ax” *ldxr 
0 
and for s > 0, define, 
If Gl cf” (ee —1+(ik- s)x) ax * de. 
0 


By the Dominated Convergence Theorem (Theorem D.3), the boundary conditions 


vanish. Integration by parts yields the following: 


J,(a@) = C(-ik — s) I Ce — 1) COs 


= ge" es (a - 1) (a — 1)a7 D1 dr, 
a—l Jo 


Another integration by parts calculation yields, 


J,(a) = os Cae = 1) (=a) + (-ik = 5) | * (ik 8)ry (0-1) gp 


The boundary terms vanish since e(~*-9)}” —1 = O(x) as  — 0. Using the 


characteristic function of a gamma probability density function: 


Comer iN. 
tka a—1,—by eee 
| T(a)” ee ( >) 
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we have, 


—1 ge S 
or, 
ik+s ave. CI(2-a) 
. = T(2— k\@ oo | a 
J(a) Cea (2 — a)(s + ik) aes (s + ik) 

Using Theorem D.3, as we let s — 0, we have, 

CI(2 — 

PBI a ae 


a-—l 


For the second integral, by using Taylor series expansion, note that for all 7,k € R 


we have 
: kr)? 
e*® _1 + ika| < ee) 
2! 
Thus, 
Cl/a ; k2 Cl/a 
i (ei —1+ ike) Caxr~*1da| < = | Cax~*"'dxr 
0 2 Jo 
Or, 
CUE yo Ke eo a 
tke : -a-1 One 5 ee a 2 : 
ii (c 1+ ike) Cax dx) < 704 |5— 55a” O(k*) 
Thus, 
T(2 
fk) =1- 0° ik +0 C=) in) + O(8) 


In the derivation of the fractional diffusion model, we will follow a similar method 
used to prove the traditional case. Suppose that {Y,,} denote a sequence of Pareto 


random variables with mean zero. These random variables can be constructed by 


letting Y,, be independent and identically distributed with X — E[X]. 


Le 


—1 
Using Proposition 1.1 for C = er and the Taylor series expansion for e*, 
—a 


as k — 0, we have, 


Ble eX -EIXD) — [1 — skELX] + (ik)* + O(k?)] |1 + tkE[X] + —(ikE[X])? + O(K?) 
= 1+ (ik)* + O(k’). 


Let S, = Y; +---+ Y, be a random walk denoting the position of the particle 
after n steps and let f(x) be the probability density function of Y,. Recall that n 
can be reinterpreted by considering the position of a particle at time t > 0 and some 
scaling constant c (i.e. let n = ct). Let u(x,t) be the probability density function 
of a random particle at time t > 0 and distance x. Since u is a probability density 


function, we require the scaling constant c~'/“ to normalize the sum S,,. Therefore, 
aie v7 ay n 
Ble the Sn) (: 2 UNE +0(%)) 
c 


Letting c > oo and n > oo yields, 


et (th) Ble **7] = ai(k, t) 


where Z is a normalized random variable. Again, it is clear that a(k,t) solves the 


following differential equation: 


dii(k, t) 
dt 


= (ik)*a(k,t). 


We refer to Example 1.1 to construct a definition for the fractional derivative. 

d® f(x) 
ro 

a function whose Fourier transform is (k)°f(k). Inverting the above differential 


For the purpose of consistency, we define the fractional derivative to be 


equation yields the fractional diffusion model provided below. 


Ou(2,t)  O°u(z,t) 
Ot  —s- «Ou 


Specific details on the inversion are supplemented in Appendix D. 


Consequentially, note that by Theorem D.1, since 
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Blethen) + eli®)” = Bfe-i*2] = fe ‘Mule, t)de 


we have, 


Yost etasey., 
oWeg, = tng 
n a 


when Z is a—stable motion with mean zero. The convergence of the distribution 
of random variables is renowned as the extended central limit theorem. 

The above fractional diffusion equation models super-diffusion. In the traditional 
case, the diffusion equation models the concentration of particles from one unit to 
another. Particles contained in a unit of higher concentration have the tendency to 
disperse into adjacent units of lower concentration. This physical behavior is known 
as standard diffusion. 

However, the fractional diffusion equation models the concentration of particles 
from one unit to many adjacent units. The concentration of particles is spread par- 
tially among all adjacent units of lower concentration. The amount of particles en- 
tering into these units are distributed exponentially. Specifically, if the unit is closer 
to the source, then there is an increased amount of particles that will enter the unit. 
This distribution of particles is known as super-diffusion. 

Compared to traditional diffusion, super-diffusion is a generalized measure of gra- 
dient change. This measure is also more realistic in many scientific applications. For 
these purposes, we tackle the inconsistencies of the Black-Scholes model by consider- 
ing its fractional counterpart. 

The solution of the fractional diffusion model is positively skewed with a heavy 
power-law tail. Specifically, as 2 + oo, then for some constant A = A(C,t,a) > 0, we 
have u(x,t) = Arv~*"' + 0(x~*"'). [7] Note that the tail behavior does not disappear 
as u(x,t) is passed to limits. This is similar to how implied volatility behaves under 
changes of maturity dates for S&P 500 options pricing after the stock market crash of 


1987. As discussed previously, as maturity increases, there exists leptokurtosis in the 
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return distribution and therefore the tails of this distribution appears fat compared 
to the tails of the Gaussian. 

Recall that in the necessary assumptions of the Black-Scholes model, we require 
a generalization of Brownian motion to compensate for the recent behavior of option 
prices. We consider a specific class of stochastic processes called Levy a—stable 
motion and readjust our assumptions so that asset prices satisfy a Geometric Levy 
motion. Given that the parameter a € (1, 2], fractional derivatives will arise and 
become crucial in the development of the Finite Moment Log Stable (FMLS) model. 
This will require some background information on a—stable motion and fractional 
derivatives. In the next section, we will first discuss about fractional derivatives and 


its interesting properties. 


1.4 FRACTIONAL DERIVATIVES 


As mentioned in the previous section, we have defined the fractional derivative as 


d® f(x) 


Be whose Fourier transform is (ik)°f(k). We will now provide an 


a function 
alternate definition of fractional derivatives using standard calculus intuition. 


Recall that the first derivative is defined by the following limit: 


U(x). fa) = fle) 


dx ho90 h 
We will assume that this limit exists. For higher order derivatives, we have, 


a f(z) _ 5 A*F@) 


dx” h>0 hn 
where A’ f(x) = )> () (—1)/ f(a — jh) for n € N. Extending this definition of 
j=0 


the difference operator for a € R yields the following: 


mM —— 
dx h-0 he h>0 he (= 


af(z) _ i; x? f(a) i Li (° 


0 
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This is known as the Grunwald-Letnikov finite difference form for the frac- 


tional derivative. We use the convention that 


(") 2 OD) 
j} j(a—j4+1) 


We will now state a proposition that shows that the two definitions of the fractional 


derivative agree with each other. 


Proposition 1.2. Let f(x) be a bounded function such that f and its derivatives 


up to ordern > 1+ a exist and are absolutely integrable for alla € R. Then, the 


Grunwald-Letnikov fractional derivative exists and as h —+ 0, the Fourier transform 

AM Fe 

op SHO) 
h= 


approaches (ik)“f(k). 


Proof. Let w; := (") (—1)/. Then, by the Binomial formula, since, 
J 


yw ed (“)ay =(1+(-1))° =0 


j=0 j=0 \J 
we have, 
Swill ( ‘yay < 00. 
j=0 j=01\I 


With the assumption that f is bounded, then for —coo < x < ov, we have the 
following uniform convergence: 
ary(e) =: (*)(-1)'s(@= 7h) <0. 
j=0 \J 
And therefore the Grunwald-Letnikov fractional derivative exists. 
We will now propose a lemma that will assist us in proving the second half of the 


theorem: 


Lemma 1.2. /f f(x) and all of its derivatives up to order n exist and are absolutely 


integrable, then for allk € R, there exists constant C > 0 such that, 


A C 
FO) < ae 


Proof. Note that for |k| < 1, since 1 + |k|" < 2, 


(1+ [AMF] <2] fe f(@)da| <2 f If(@)lde. 
Therefore, 
f(A) < ape | eae. 
Alternatively, note that by Example 1.1, we have, 
Fm) = Ry fe f(a)ae, 
For |k| > 1, since 1 + |k|” < 2|k|”, 


(1+ |kI")1F(&)] < 21k|" 


Go le” FO (a)da| <2 [|f(a)lar. 
Therefore, 


P 2 fy em 
\f(k)| < eerie ii (o) \de. 


Letting C := max {2 fede. 2 f f(x) ae} proves the lemma. 
Using Lemma 1.2, it follows that for all k, 


|(ik)* F(k)| < 


A 


Given that n > 1+ a, (ik)*f(k) is absolutely integrable. By applying Theorem 
D.2, we can deduce that there exists a function with Fourier transform (ik)“f(k). By 
an earlier discussion, we have denoted this function as the fractional derivative. 


Since, 


fe Fe —alde = fe F(y)ay =e fe F(y)dy == ef (h), 


Ze 


the Fourier Transform of A” f(x) is 


foot ne Be fone 


j=0 


3 (“Jn Jeti" F(A) 


j=0 \J 


= (1—e"™")*f(k) 


And the Fourier Transform of h~*A" f(z) is: 
See ee _ e 7 tkh OF 
nad =e 2F(a) = arecinny® (7S) FG) 
= (ik)* (1+ O(ihk))* f(k). 
Letting h > 0 yields that for all h, 


/ etn A" F(a)dax > (ik)*f(k). 
el ) d* f(x) 


Since the Fourier transform of converges pointwise to 7 , the conti- 
po 


nuity theorem for Fourier ee. proves the theorem. 


Remark: The term w,; defined in the previous proof is known as a Grunwald 


weight. Note that by definition and properties of I'(-), 


os (F\a =» Aes. — =e ae) 
i Tg+l(a-j+1) TYyt+i1rd-a) 


Applying Stirling’s formula P(2 + 1) ~ V2rx (=) as x — oo yields, 
e 


: r j—a-1 
tires inant (heat) pra leetl, 
7 T=a)Vo 5 j 


Letting 7 — oo yields, 


Ta -—a—1 
“a rd - a)” 
Note that, 
Aa e0 
Ree = Aa)" | Fle) +X wy fle— 580 
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Let 0 < a < 1 and for j > 1, define b; = —w; so that, 


bye pqaayh 1 I 00 and 8) = 1 
Then, 
a = (ag? Li y bj + S(-8))F(e — jAz) 
= (Ar) SUsla) — f(a — jAz))bj 
© DR) — Fle jae) yay ae 
= [Ue fe- Mpg grew 


Passing to limits yields, 


d® f(x) 


dx@ 


= [if@) - fle- Wp 


After applying integration by parts, we have 


de f(z) 1 oo ss 
dz faa dat * — WY dy. 


Rename u = x — y and taking the derivative outside yields, 


df (x) 1 od i 


f(u)(a — u) “du. 


ae T(1 — a) dz J-c 

The above form is renowned as the Riemann-Liouville fractional derivative 
for 0 < a < 1. Analogous forms can be constructed for other values of a. This 
definition is necessary for the derivation of the Fractional Black-Scholes equation sat- 


isfying the FMLS model. We will eventually restrict our attention for all a satisfying 


l<a< 2. 
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1.5 a-STABLE VARIABLES AND FMLS MODEL 


In this section, we will develop a framework for a—stable random variables and discuss 
about some properties they hold. We will further use these properties to construct the 
Finite Moment Log Stable model. Readers can refer to Samorodinitsky and Taqqu 
for proof of these properties. [12] 

Firstly, we propose the definition of a Levy process to construct a similar definition 


for a—stable motion. 


Definition 1.2. A time-dependent random variable X is said to be a Levy process 
if and only if X has independent and stationary increments. The Levy-Khintchine 


representation of the log-characteristic function of X its provided below: 


1 
log Efe**)] = mitk — <0?tk? +t the 1 _ ikadiy<)dW 
og Ele |= 57 + aoa tkaljz\<1) 


with driftm € R, volatility o > 0, indicator function I and Levy measure W 


satisfying: 
: min{1,27}dW < oo. 
R\{0} 
If Levy measure dW is of the form w(x)dx, we define w(x) to be the Levy density. 
For our purposes, the Levy density of a—stable process is defined as: 


bs 1(1—g)D|2|\--*@ a <0 
WLS\t) = 
$(1+ 8)Dx1* oo 


where D > 0, and skewness —1 < 6 <1, and0 <a < 2. Applying Definition 1.2 


for this Levy density yields the following characteristic exponent V(k): 
; deter ae i TO 
W(k) = ik — 5 lf o (1 — i8(sgn(k)) tan =) 


with drift w and volatility o. This motivates us to study a class of Levy processes 


called a—stable processes. We will shortly show that this class of stochastic processes 


25 


satisfy the desired conditions reflected in the S&P options market. We provide the 


following definition: 


Definition 1.3. A random variable X is said to have an a—stable distribution if for 
the given parameters, we have the following conditions: index of stability1 <a <2, 
drift u, dispersiono > 0, and skewness parameter —1 < 6 < 1, then the characteristic 


function of X, yx = px(k), ts: 


ae Bet] = exp ik — |k|%o* (1 — iB(sgn(k)) tan =)| 


Letting a = 2 yields the characteristic function of a Gaussian random variable W 


with mean pz and variance 207. Specifically, we have, 


~w = 


E[e*”] = exp liky — k’o?|. 


In this class of stochastic processes, the thickness of the tails remains invariant 


under time. This is shown with the following property. 


Property 1.1. Let X be a—stable with a < 2, dispersion o, skewness 3, and drift 


bu. The tails remain "fat" and the tail probabilities are given as follows: 


1+ 
Jim ASP (AEX. A) Cla) 5 Bae 
where, 
| ee 
sau -1 oF 1 
C(a) = (/ ae) _ ) T(2—a)cos(ra/2) 
2/0 a=1 


When X is maximally negatively skewed with zero drift (i.e. B = —1, 1 =0), we 


have 


; y  \ral2a=2) y \alle=n) 
page 27a(a — 1) aa) re [to tl (acixa)) 


where, 
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~1/o 
Clo;a) =a (co eo) 

This property also suggests that for a < 2, the tail probabilities behave like \~° 
and therefore exhibits a "power-law"-like behavior. Note that since the tail behavior 
does not flatten out as maturity increases, a—stable variables violate the implications 
of the central limit theorem. This suggests that variance for log return is infinite. 
The questions of the existence of a martingale measure and whether option values 
are finite or infinite arise. 

By setting the condition of a—stable motion to have maximum negative skewness, 
we claim that conditional moments of all orders exist for 1 < a@ < 2. To show this, 
we need to propose three additional properties of a—stable motion that will assist us 


in the proof of our claim. We will first start off by providing a definition and some 


notation. 


Definition 1.4. Let Y be a random variable with probability density function f(x). 


We define the Laplace transform of f(x) to be the following: 


f(s) = B{f(@} := Ele-*] = i * ent f(t) dt. 


We define the two-sided Laplace transform of f(x) to be the following: 


f(s) = B{f(t)} := Ele-*] = 3 en" f (t)dt. 


Similarly to how we have notated W, = W(t) for Brownian motion, we also 
propose a notation for Levy a—stable motion. Formally, we will denote L@? = 
L(t) as the standardized Levy a—stable motion with tail index 0 < a < 2, and 
skew parameter —1 < @ < 1. Furthermore, by Definition 1.3, given that there are 
four parameters, we will conveniently notate the distribution of a—stable motion as 
Lalu, 7,8). We say that if X has a stable distribution with the above parameters, 
we notate this as X ~ La(p,¢, 2). 


We will now proceed to mention some additional properties: 


yas 


Property 1.2. Let X ~ L,(",0,8). The two-sided Laplace transform of X is not 


finite unless 6 = 1. For 8 = 1, we have the following Laplace transform: 


E[e**] = exp (—sn — s%o% sec *). 


Property 1.3. For any 0 <a < 2, 


X ~ L{0,0,8) 4 —X ~ L, (0,0, -8). 


Property 1.4. For any0 <a < 2, we have E| 


| 


X|?| < co for any 0 < p< a and 


X|?| = co forp>a. 


What makes a—stable processes so appealing for the construction of our new 
model is described by these properties. As mentioned previously, Property 1.1 sug- 
gests that the tail behavior is invariant to maturity. We say that this property 
exhibits self-similarity. Property 1.2 defines the Laplace transform for a—stable 
random variable X with maximum positive skewness, and Property 1.3 allows us to 
determine a closed form for the Laplace Transform by relating X to —X. Lastly, 
Property 1.4 describes the finiteness condition for the moments of X depending on 
Q. 

Let L°’ be a Levy a—stable process with skew parameter 3. Define S; to be the 
stock price and assume that 5S; satisfies the following stochastic differential equation: 


for time 0 < t < T, index of stability 1 < a@ < 2, and volatility 0 > 0, 


Bt. (r — q)dt + odLo" 
St 


where r and q respectively denote deterministic parameters corresponding to the 
risk free rate and dividend yield. We selectively restrict 6 = —1 to obtain finite 
moments of S; and negative skewness in the return density distribution. Specifically 


by Properties 1.3 and 1.2, for n > 0, 


[exp (noL?')] = exp (-tn0" sec _ < 00. 
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This suggests that a martingale measure exists for 6 = —1. Furthermore, we 
restrict a to satisfy 1 < a < 2 so that S; remains unbounded. We define the 
above stochastic differential equation to be the Finite Moment Log Stable Model 
(FMLS model). In regards to this model, we state and prove the following proposi- 


tion. The proposition is two-fold. 


S 
Proposition 1.3. Let s, := log =F be the log return over maturity 7 = T —t. Then, 


St 
TO 
i) 8, ~ Lo((r—q4+ p)t,07!/%, -1) with convexity adjustment ps = 0% sec me 


ti) For alln > 0, the nth conditional moment of S'p is well defined and given as, 


Bc 97 SP (nr — q+ p)T — T(no)* sec *) 266 


Proof. Let S; satisfy the above stochastic differential equation for all 0 < t < T, 


1<a<2,andoa>0. We can re-express S'r in the following exponential form: 


Sr = Spe"-97 exp (por + oLe1) 


where u is specifically chosen so that El (qr + gig) = 1. Note that by Property 


1.3 and 1.2, we have 


Elexp(oL¢~")| = Elexp(—oL?")] = exp (-r0" sec *). 


This requires 4 = 0° sec > < oo. With this choice of uw, we have the following 


exponential form, 
Sp = Sye-T#)" exp (9 L&—?), 


Therefore, the log return s, satisfies 


S: 
5 = log “A =(r+qtprt+oL*. 
t 


Note that s, is a—stable distributed with mean (r—q+y)7, dispersion or/and 


skewness 6 = —1. We can further deduce that by Property 1.4, the variance of s, or 
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any moments of order higher than a is not finite. This proves the first statement of 
the proposition. 


Note that for n > 0, by Property 1.2, 


E[S2] = Serato) lexp (no Le") ] 


= S? exp (n(r —q+p)t —T(no)* sec =) < 0O. 


This proves the second statement of the proposition. The proof is complete. 


With the proposition proven, the FMLS model performs better than the standard 
Black-Scholes model in reflecting how S&P option prices behave after the 1987 U.S. 
stock market crash. As desired, under the FMLS model, the return distribution 
contains "fat" tails with the extra condition of infinite return moments and finite 


price moments. 


1.6 DERIVATION OF BLACK-SCHOLES EQUATION UNDER FMLS MODEL 


In Section 1.2, we have derived the Black-Scholes model and its corresponding par- 
tial differential equation using the assumption that asset prices satisfy a geometric 
Brownian motion. Alternatively, by applying Ito’s Chain Rule (Theorem B.1), the 
Black-Scholes model states that for asset price S;, the log price satisfies the following 


stochastic differential equation: 
1, 
dibes\= ¢ -q-50 ) dt + odW, 


with risk-free rate r > 0, dividend yield g > 0 and volatility 0 > 0. The price 
of a European call option C(S,t) with S; satisfying the above stochastic differential 
equation is given in Section 1.2. With a change of variables 7; = log S;, the Black- 
Scholes equation can be rewritten as the following advection-diffusion type equation: 


(x - 50°) gOS =1OGi st) 


OC (et) lk sOC(Et) 
+ ~O 
Ox 


ot 2 Ox? 
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with C(x,t) representing the price of a European call option under the Black- 
Scholes model and risk-free rate r. 

In this section, we generalize the derivation of the partial fractional differential 
equation that arise from the FMLS model. More details can be referred in Cartea 
and del-Castillo-Negrete. [4] We readjust the Black-Scholes model by assuming that 
the stock price follows a geometric Levy process with maximum negative skewness. 


This motivates us to consider the following stochastic differential: 
d(log S;) = (r —q—v) dt +o0dL?"'. 
The above equation has the following solution: 


T 
Sp = Stexp ( —q-—v)t+ o | it) 


where 7 = J —¢ and v is a convexity adjustment so that E[S7] = e975}. 


In determining the corresponding fractional partial differential equation that models 
the behavior of the value of a European call option V(z,t), we prove the following 
proposition that provides a differential equation which governs the behavior of its 


Fourier transform. 


Proposition 1.4. Let V(x,t) be the value of a European call option under the con- 
vention that x := x(t) = log S;, for asset price S; withO <t < T. We assume that 


S; satisfies a Geometric Levy distribution. Then, the Fourier transform of V(z,t), 


A 


V(a,t), satisfies the following differential equation: 


av (k,t) 


as (r +ik(r —v) — U(—k)) V (k, t) 


where V(k) is the characteristic exponent for some Levy process. 


Remark: For the FMLS process, by Definition 1.2, we have, 


vO=2Gre slkl%o* (1 shot ie *). 


dl 


Equivalently, we have, 


U(k) = —i0* sec > (l= 6) Gk)? + (1+ Bye tk)*). 


Thus, for @ = —1, we have, 


1 
W(—k) = a sec > (-ik)*. 


More details are provided in Benson, Meerschaert, and Wheatcraft. [1] 


Proof. Assuming the market is risk-free, the value of the European call option can 
be written in terms of the expected value of the final payoff. Let II(x,7’) be the final 
payoff of the portfolio. Then for interest rate r and maturity 7 := T(t) = T — t, we 


have, 


dV (a,,t) = rE[I (a7, T)|dt = —rE|I (x7, T)|dr. 


Equivalently, 


V (2, t) = e "EM (a7, T)]. 


Assuming that II(2, 7’) has a Fourier transform, applying Theorem D.2 on II (x7, T) 
yields, 
—r(T-t) 
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e 


Vege) = 


E | / * enter fl(k, T)dk. 


By linearity of expectation and using the characteristic function of xr = log Sr, 


we have, 


V(x,t) = © ie i . E (e**") T(k, T)dk 


eT (Tt) oo ; ; eae 
= / e ika,—ik(r—v)(T t) (7 t)w( NTI(k, T) dk. 


Applying Theorem D.2 on V(x,t) and equating forms yields, 


V (k,t) = r—ik(r—v)+W(—k)|(T OT1(k, T) 
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with V(k,T) = II(k,T). It can easily be checked that the above function solves 


the differential equation: 


aV (k, t) 


ae (r + ik(r — v) — U(—k)) V(k, 2). 


By the previous remark, we have established that under the FMLS model with 


maximum negative skewness, 


1 at 
W(—k) = —x=0% sec —(-ik)*. 
(—k) 57 SECS (—ik) 
; ; : Lat SOR oo oe 
With convexity adjustment v = = 50 Sec a in Proposition 1.3, by Proposition 


1.4, we have, 


DUCA ae Oe a Re Lea 
oy = (r+ ik (r+ 50 sec =) iat SeC mee) ) (6,0). 


Using Theorem D.2 on both sides yields, 


OV(z,t) | Niet OR VOM (GE) De OL V GSE) « 
BE (x t 5 sec ) oe a ae = TV (0.6) 
where, 
O°V (a, t) _ 1 d m -—a 
dre —s«dT(1—a) dz f oe 


The above fractional partial differential equation is renowned as the Black-Scholes 
Equation under the FMLS model. This is the equation that we will be studying the 
solutions for throughout this paper. 

Remark: Note that if ~ = 2, the derived Black-Scholes equation degenerates into 


the classical Black-Scholes equation in its advection-dispersion form, given below: 


OC(z,t) | 1 20°C(z,t) : (r- 5) OC (a, t) 


Ox 


Ot 5 72 =a, b): 
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CHAPTER 2 


ANALYTIC SOLUTION 


From the previous chapter, we have established that under the FMLS model, the 
log price under the risk-neutral measure satisfies the following stochastic differential 


equation: 
d(log S;) = (r —q —v)dt + odL?' 


where r and q are the risk-free interest rate and the dividend yield respectively, 
and v = —50° sec > represents the convexity adjustment. Let V(z,t) be the price 
of the European call option with « = 2; := log S;. From the previous chapter, we have 
shown that V(x, t) satisfies the Black-Scholes equation under FMLS given below: for 


all. <P andl <or 2. 


OV (a, t) 1 am\ OV(az,t) 1 at O°V (a, t) 
t+ (r+ -o* — 9% sec“ = rV(a,t 
ey ( on" see > ) ae 57 sec Ana rViest) 
with boundary conditions, 
max{e” — k,0} for European call option 


Vast) =H 
max{K — e”, 0} for European put option 


where K is the strike price and, 


OPV (at) 1 d ‘: 
Ox (1 — a) dz 


V(u, t)(a — u) “du. 


In this chapter, we will determine a closed-form analytical solution of the Black- 
Scholes equation using Fourier transforms. We will state a few theorems relating our 
solution to the Black-Scholes pricing formula provided in the previous chapter. More 


information can be referred to Chen, Xu, and Zhu. [5] 
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2.1 PRELIMINARIES 


In our derivation of the generalization of the Black-Scholes formula, we require addi- 
tional definitions and properties of a—stable motion. We also require a lemma that 
will help us prove a theorem discussed later in this thesis. 

By Definition 1.3, the characteristic function of a—stable random variable X with 


drift 4, dispersion o > 0, and skewness parameter 3 € [—1, 1] is 


p(k) = exp (ihn — |k|%o° (1 — iBsgn(k) tan =)). 


Letting 6 = 0 yields a symmetric distribution with translation constant w and 
scaling factor o®. This can be shown by considering the Levy density function. 


Namely, for 6 = 0 and for some constant C' > 0, the Levy density becomes 


Claire: -grenG 
wis (2) = 
Ca zt >0 


Note that the Levy density function under this choice of @ is symmetrical. Since 
wrgs correlates to the behavior of the distribution, it implies that the Levy distribution 
is centered and symmetrical. Therefore, for 1 < a < 2, |y(k)| = exp (—|k|®). 


In general, under new centring constant 6 satisfying 


| < a 0O<a<l 


2-a l<a<2 


and by removing the drift and dispersion by letting 4 = 0 and a = 1, the log 


characteristic function of X, W(k), is defined as follows: 


W(k) = —|k|° exp (:Fsn(e). 


This motivates the following property. 


Property 2.1. Let fog(x) be the probability density function of a—stable variable 


X. fag is the Fourier transform of characteristic function ~p(z), given as follows 
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fa,a(X) = ~Re [~ exp (ive =~ 2° exp (=) az. 


Now, we propose the definition of a function that plays a central role in the 


determination of solutions for fractional partial differential equations. 


Definition 2.1. Let A;,B; > 0,0 <m<qand0 <n <p be constants. Let a;,b; 
be complerz numbers such that no pole of 1(b; — Bjs) for j =1,--- ,m coincides with 
any pole of ((1—a,; + A;s) for j =1,--- ,n. We define the Fox H-function to be 


the following function: 


m,n min ee ( 
yg (2) = Aye lz 


~ Oni =e 7 


j=n+1 Tenn 


wo je; Ass) se 


Ie 
Tf paigvh k=); +2878) 


where C' is a contour in the complex plane such that and “= lie to the 


atk 
B; 
right and left of C, respectively. 


The probability density function of stable variables can be rewritten in terms of 
this special function. We will propose the analytic form of the probability density 


function using the Fox H-function. 


Property 2.2. Fora > 1, the analytic form of fa, is given as follows: 
_ eae 
fo,e(%) = me av 


Soon, we will see that the closed-form solution of the Black-Scholes equation can 
also be written in terms of the Fox H-function. Crucial to the understanding of Fox-H 


functions, we will propose the definition of Mellin transforms: 


Definition 2.2. We define the Mellin transform of a function f to be the follow- 


ing: 
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M{f(a)} = [oo f(a)de. 


The inverse Mellin transform is defined as: 


MLE) = = fo 2s) 


where the line integral is taken along a vertical line in the complex plane. 


Lastly, we propose the following lemma that emphasizes the properties of the 


Gamma, Function: 


Lemma 2.1. (Gauss Multiplication Formula) For all z € {-= :m¢ Nn}, 
n 


n—1 
I[v (: + =) = (2r)"2 n2-"T (nz). 
k=0 Ms 


Proof. The Euler form of the Gamma Function is provided below: For all z € R\Z, 


ae I (+z) (+2) = be oy GE 


n=1 


Using this definition of the Euler form, Stirling’s formula and properties of the 


Gamma function, namely ['(z + 1) = zP(z), we have: 


eee eee) 


mlmZztk/n-1 


I 


Van (m2) metk/n—32 
— li : . 
modo (nz+k)(nz+k+n)---(nz+k—n+mn) 


Therefore, making a substitution of mn +> m and reapplying the Euler form of 


the Gamma Function, Stirling’s formula, and properties of the Gamma Function, we 
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have: 


n moo (nz) (nz+1)---(nz-—1+mn) 
(Vm) (2) 
moo (nz) (nz +1)---(nz -—1+ mn) 
(v2r)" (2) mm Vn t/2—ne 
moo (nz) (nz + ) -++(nz—1+m) 


(v2n)" milmn2z—lyl/2-nz 


moo (nz) (nz+1)---(nz-—1+™m) 


— 
ee 


= (Qn) PVP pV2-m2P (yz) . 


The proof is complete. 


2.2 CLOSED-FORM ANALYTICAL SOLUTION 


1 
Firstly, define rT := ae sec SL —t). The Black-Scholes equation becomes the 


following: 


with boundary conditions V(z,0) = II(a) for H defined previously and y = 
ant ; a ; 

Tz O08 SY can be reinterpreted as the relative interest rate with respect to 

volatility of fractional order. Taking the Fourier Transform on the above fractional 


partial differential equation yields the following: 


with boundary conditions V(k,0) = II(k) where II is the Fourier Transform of the 


payoff function I]. The above differential equation has the following solution: 


V(k,r) = e-VT(k) exp (—(1 — y)rik — |k|“r). 
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Applying the convolution theorem for Fourier Transforms yields the following 


convolution, 
V(z,T) =e 7V (a, 0) « F-‘{exp (—(1 — y)rik — |k|%r)} 


Define P(x) := F-'{e7""7\. Applying the shift theorem for Fourier transforms 


yields the following: 
V(a,T) =e 7 V(az,0) * P(a — (1 —)r7) 


By previous discussion, e~!*'" is the characteristic function of a centered and sym- 
metric Levy distribution. Therefore Property 2.1 implies that the Fourier inversion 
of e~!*I"7 is a multiple of the Levy stable density function fo. Applying Property 


2.2 yields the following representation of P: 


ws qi/a Fou Tia} gla" ?? | r/o ip a 
(0, 1) ee s) 


Thus, 


aes) qzi/a qci/a 


Vea) = ye eT Kk) : foo ( moa oe vel) dk 


Or equivalently, 


V(2,7) = ‘| * eM (2 — (1 — 9) — 7°) fxg (|) ab. 


—oco 


Let V,(xz,7) and V.(%,7) be the price for a European put option and a European 


call option respectively. For European put options, we have: 


V,(a,7) = Ke foo (|kl) dk — ay exp (—7 — r/*k) fo ({hl) dk 


—log kK —-(1-— 
where d, = eee 7 ( ae An analogous solution for European call op- 
qri/a 


tions can be determined similarly. We will refer the above formula as the Black- 
Scholes Formula under the FMLS model. In the next section, we will discuss about 


the properties associated with this formula. 
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2.3. PROPERTIES OF BLACK-SCHOLES FORMULA UNDER FMLS MODEL 


Recall that when a@ = 2, the FMLS model degenerates to the traditional BS model. 
Furthermore, for boundary conditions 0 < t < T and 0 < S < co such that C(S,t) ~ 
Sas S > oo, and C(S,T) = max{S — K,0} for strike price AK’, the solution of the 


Black-Scholes equation is provided below: 


C(S,t) = SN(d,) — Ke") N(dy) 


where, 
N(o)3= ae exp (~5s*) ds 
and, 
1, _ 8(S/K) + (r+ 38°(0 = 8) 
: Birt 
leg (S/R) = ser 5) 
dy = B/Pot =d, — BVT —-t. 


We will show that by setting a = 2, our derived Black-Scholes formula degenerates 
into the standard Black-Scholes formula mentioned above. 

Note that the Fox H-function is defined as a Mellin transform of a rational ex- 
pression of the Gamma function. We will use this fact in the proof of the following 


theorem. 


Theorem 2.2. For European put options, the Black-Scholes formula under FMLS 


model degenerates to the Black-Scholes formula as a — 2. Equivalently, 


lim, V,(z, 7) = Ke~ N(—do) — e” N(—d,) 


x —log kK +(y-1)r 
V2r 
NG) ae exp (-53) ds. 
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, dy = d, — V2r, and 


where d, = 


Remark: The above variant is equivalent to the standard Black-Scholes formula. 


This is achieved by making the following substitutions: r := y, 6 := V2, and S := e”. 


Proof. By Property 2.2, we have, 


foo (Il) = lim fa.o (Al) 


2G Is, VE 
= lim —H2} | | ( aie) (oe) 
(1) x5) 
idl 1d 
= -H}} |m| (Sam) ma 
(OpL)e Gs) 
if 
1 999 
= Lan || (2? 
(0, 1) 
where, 
a1,A if ['(b) + B 
isEee) -= Hy? |z (a1, Ar) = _f a PS) sds, 
(b1, By) mi Jo D(a, + Ais) 


The Fox H-function is defined as a Mellin transform of a rational expression of 
1 [(s) 


the Gamma function. By Definition 2.2, we see that M{ foo ({k|)} = ==4-—. 
Applying Lemma 2.1 for n = 2 yields: 
I 1 loz 
[(2)P (- i 7 — (2n)223--T (22) 
or equivalently, 


P(2)P (< zs 5) = 91-2 FED (22), 


Substituting z= ; yields, 


Therefore, 


Taking the inverse Mellin transform yields, 
1 I(s ge eK /4 
fro (el) = at | a at | a) 
2G +38) 4,/m 2/m 


This shows that fo (|k|) is identical to the Gaussian density function. Thus, by 


the Black-Scholes formula under FMLS model, with some algebraic manipulation, 


the Black-Scholes formula is obtained. Thus, the theorem is proved. 


We can also determine the asymptotic behavior of our newly derived solution of 
the fractional partial differential equation. This will allow us to study the behavior 
of option pricing for extreme values. We will now propose the following theorem and 


provide proof. 
Theorem 2.3. For European put options, for 1<a< 2, we have: 
Se _ 
_jim. V,(2,7) = Ke” and lim V,(x,7) = 0. 
Proof. Note that as 7 — —oo, we have d, — —oo. For symmetric Levy density, we 
have / fo. (|k|) dk = 1. Thus, we have 
lim Vp(e.7) = Ke” [foo (lhl) dk — tim. e* [exp (7 — 7/*8) foo (Ih) dk 
mie": 
This proves the first statement of the theorem. To observe the behavior of the 


Black-Scholes equation under FMLS model for European puts as log price approaches 


infinity, we note that 


lim V,(e,7) = Ke-? [ fuo(|kl) dk — lim e [exp (—r —7/*k) fa (\RI) dk 
zoo «=P 8S 2 LOO dy : 


== lim : exp (-r — leg) foo (\k|) dk. 


roo E-# 
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Applying L’Hopital’s Rule and using the fact that any density function approach 


ZeTO as  —> Oo, we have, 
jim V,(z,7) = A lim: ete Me Fo (|dil) 
= Ke tim fo (Id) 


=0 


since d; + oo as x > oo. This completes the proof of the theorem. 


The implications of Theorem 2.3 are expected. Under the FMLS model, a Euro- 
pean put option is expected to hold the current value of the discounted strike price 
if asset price becomes really small. On the other hand, a European put option is 
expected to hold no value if the asset price increases to infinity. 

In our derivation of the traditional Black-Scholes model, we specifically derived 
a partial differential equation that determines the behavior of European call options 
over maturity. Our derivation of the Black-Scholes equation applies to the behavior 
of European put options. It is our interest to determine if the relationship that lies 
between put option pricing and call option pricing satisfies the put-call parity. We 


will state and provide a proof of the following theorem. 


Theorem 2.4. For 1 < a < 2, under the same strike price K and maturity T, the 
price of a European call option and the price of a European put option satisfies the 


put-call parity. Equivalently, 
V.(2,7) — Vp(z, 7) = e® — Ke”. 


Proof. Given that the price of both European call option and put option satisfies the 
Black-Scholes equation under FMLS model, for V._,(#,7) := V.(a,7) — V,(#,7), we 
have, 


OVers 
OT 


OVS 
Ox 


Or Vag tet) 
Ox® 


=(7-1) — YVe-p = 0 
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2 
with boundary conditions V._,(#,0) = e” — K and y = aoe cos Using the 
o 
same approach highlighted in our derivation of the Black-Scholes formula, we have 
Vep(t,t)=—Ke-™ | fao (lhl) dk +e" [exp (= —7/*h) faa (kl) db 


In the proof of Theorem 2.3, we have established that / foo (\k|) dk = 1. 
Furthermore, upon realizing that the second integral is the Fourier transform of 


foo = Fe" 1}, we have: 


—oo 


e - exp (-r — qk) fanliehdhseo* a exp (-i (-ir'/*) k) fo (\k|) dk 
= ert (-ir'/*) 


a—T+[i(—ir!/%)]¢ x 


=e = €°. 


Thus, we have, 
Veni = he. 


And the theorem is proven. 


In this paper, we have proposed the Black-Scholes formula to determines the price 
of a European call option. To provide consistency of this paper, we will conclude this 
chapter with the generalized analytic form of the Black-Scholes formula for European 
call options by applying the previous theorem. 

The Black-Scholes formula for European call options is provided as follows: for 
log price of a European call option « and time T = —50” sec ee Gy — t) satisfying 
0<t< TT, the price of the European call option is analytically given as follows: 


dy 


dy 
Vew,7) =e* f° exp (-7 = 7k) foo (Ih) ak — Ke [fag (\kl) dk 


—log K —-(1- 2 
8 ; eos and y = cos and 
i/o oe 


where d,; = 5 
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“—— 
SS 


1 
Ha’ \|k| 


fai (|A\) ae 
> (0, 1) 


La 
NIE Ne 
NleF Nle 
Na 


when AS5 is the Fox H-function defined in this section. 


Remark: Note that by put-call parity, we have, 


ic exp (-r - rok fa, (kl) dk = 1. 


—co 
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CHAPTER 3 
NUMERICAL METHOD 


In this chapter, we will use the Finite Difference Method (FDM) to determine solu- 
tions for the Black-Scholes Equation under the FMLS model. We will also propose 
theorems that highlight the stability and the convergence of the method. 

Let V(a,t) be the price of a European vanilla option with respect to log price 
x and time t. Let r be the risk-free rate of interest, 0 > 0 be the volatility, and 
1 <a < 2 be the fractional parameter. Recall that the Black-Scholes Equation is 


given as follows: 


OV(z,t) | ha, OT OV (Et) 1 4. om OPV (at). 
BE (x + —o sec ~) a 57 860 aw TV (at) 
where, 
OV and) 1 d st - 
= —u) “du. 
Ox T'(1— a) dx i me a 


The boundary conditions are given as 


max{e” — K,0} for European call option 
Viet) =e) 
max{ Kk — e”,0} for European put option 
where K is the strike price. 


In Chapter 2, we have applied the following change of variables: 
1 
= a sec mao —t). 


2 


This resulted in the following variant of the Black-Scholes Equation: 
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Oa 7 ioe OV aT) 
Or zs Or Ox® 


with boundary conditions V (x, 0) = II(a) for II defined above and y = -* cos = 
Throughout this chapter, we will solve this variant since this will simplify conditions 
for the finite difference method that will be specified in this chapter. 

Note that the log price of the asset x satisfies —oo < x < oo. To apply the finite 
difference method, we truncate the interval and assume that for some constant y > 0, 


the log price x satisfies -p < x < y. The constant y will depend on the maximum 


price of the asset and can be readjusted appropriately. 


3.1 STABILITY OF METHODS WITH STANDARD GRUNWALD APPROXIMATION 


In this section, we develop a scheme using standard Grunwald weights to approximate 
the derivative of fractional order and establish that by using this approximation, both 
implicit and explicit Euler methods yield unstable results. [8] 

Let h = ae x; = —yptth for i = 0,--- ,m, let T, = nAr where Az represents the 
incremental change of time in the interval 0 <7 < —50" sec ao for n = 0,--- ,m, 
and’ Vi Aas rp) 

The discrete Grunwald approximation for the fractional derivative is given by the 


formula: 


OV ey 2a 1 QT(k-a) 
T(k+1) 


V(a —kh,r). 


T(k — a) 
I(-—a)P(k + 1) 
weights. We analyze both forms of the Euler method to determine if the methods are 


be the normalized Grunwald 


To simplify notation, let g, := 


stable. 


Theorem 3.1. The explicit Euler method of the Black-Scholes Eqauation under 
FMLS model with the discrete Grunwald approximation to the fractional derivative is 


unstable. 


AT 


Proof. Under the explicit Euler method, the Black-Scholes Equation reverts to the 


following: fori =1,---,m—1,andn=1,---,m—l, 
Vera Ve VPs 2 lle 
t We fe es vr =f t i—1 Vr 
Ar VG + (y ) h pa 22 98 i-k 


Solving for V,"*" yields, 


=f Ngee 
yr _ uve = (7 aM =) ATVs + ve So ORV = VP Ar 
k=2 


A At 
= —1)+ Fe Assume that the following numerical value only 


have error, i.e. V,2 = V,° + €? for some €? > 0 dependent on i. This error produces a 


where p = 1+ 


perturbation on the numerical value of V,! = V,! + et. Then, we have 


0) = nh - (T+ Fe 


Ft pe) AaWi + FED ovine — Wer 


This implies that ¢; = pe?. Recursively, we can deduce that €” = y"e?. For the 
explicit Euler method to be stable, we require the error to be sufficiently small. Thus, 


js roust satisfy |u| <1. Note that, 


A 
1-2 = a. = |p|. 


This implies that the explicit Euler method is unstable. The theorem is proven. 


Theorem 3.2. The implicit Euler method of the Black-Scholes Eqauation under 
FMLS model with the discrete Grunwald approximation to the fractional derivative is 


unstable. 


Proof. Under the implicit Euler method, the Black-Scholes Equation reverts to the 
following: fori =1,---,m—1andn=1,---,m-—1, 


Vert — VP 
AT 


vrttioymt 4 a 
i 1 4 a D 9ViEE 


=0 


SE ea) 
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Solving for V,"*" yields, 


oe 
A MIN 
where v = (: — (4 —1)- =| . Again, assume that the following numerical 


value only have error, i.e. V,° = V,° + €? for some €? > 0 dependent on 7. This error 


produces a perturbation on the numerical value of V;! = Vj + et. Then, we have 


~ a = 1 1 a 
Vi = vVo = (-1e" = Y Ve % De aves] Ar 


SOV PG Ve Eve, 


a 


This implies that ¢} = ve?. Recursively, we can deduce that ¢” = v"e?. For the 
explicit Euler method to be stable, we require the error to be sufficiently small. Thus, 


y must satisfy |v| <1. Note that, 


(1-6-0 - 3) 


This implies that the implicit Euler method is unstable. The theorem is proven. 


1l< 


Therefore, both explicit and implicit Euler methods yield unstable results and the 
discrete Grunwald approximation cannot be used to solve the Black-Scholes equation 
under the FMLS model. In the next section, we will make a slight adjustment to the 
Grunwald approximation of the fractional derivative that will allow the Euler method 


to be consistent and unconditionally stable. 


3.2 STABILITY OF METHODS WITH SHIFTED GRUNWALD APPROXIMATION 


We will begin this section by defining the shifted Grunwald approximation of the frac- 
tional derivative. With the same conditions as before, we define the shifted Grunwald 


approximation as follows: for positive integer p, 
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1 1 AT(k-a) 


AnV = V(a —(k- : 
h (2, T) I'(-—a) he Xu [(k ate 1) (x ( p)h, T) 
Furthermore, let AV(z,7) = wv et) be the standard Grunwald approximation 
Ss pad 


of the fractional derivative. We claim the following theorem: 
Theorem 3.3. As h > 0, A,pV(x,7) = AV(a2,7T) + O(h). 
Proof. Taking the Fourier transform of the definition of the shifted Grunwald approx- 


imation, we have, 


m 


AV (7) = Fe "(Jerr ,7 


=, J ihe _ eh yey (k, T) 


he 
1 — ptkh\ oe | s 
= 5,(-ikh)® (Se eV (k, 7) 


= (—ik)*w(—ikh)V(k, 7) 


1—e*\* 
= =1- (p- ) z+ O(|z|*). Since for all z, there 


exists C' > 0 such that |w(—iz) — 1] < Cz], we have, 


with w(z) = e”? ( 


AnV (kT) = (—ik)°V (k, 7) + (—tk)*(w(—-ikh) — 1)V(k, 7) 


= AV(k,7) + o(h, k,7) 
when $(h,k,7) = (—ik)*(w(—ikh) — 1)V(k,7). This implies that 
|9(h, k,7)| = [kI°C|ak| Vk, 7) 


Applying Theorem D.2 implies that as h > 0, A,V(x,7) = AV(az,7T) + O(A). 


Theorem 3.3 discusses how as h — 0, the shifted Grunwald approximation behaves 
similarly to the standard Grunwald approximation of the fractional derivative. This 
enables us to consider finite difference methods using this altered approximation. We 


will propose a theorem stating our observations. 
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Theorem 3.4. Let 1 <a<2. Then the implicit Euler method of the Black-Scholes 
Equation under FMLS' model with the shifted Grunwald approximation, given as fol- 


lows, 


V(a — (k —1)h,7) 


2 
jor hs a to the fractional derivative is consistent and unconditionally stable. 
m 


Proof. Note that there exists log price C' such that V(C,7) = 0. Let —py < C. The 
initial condition V(—y,7) = 0 implies that for all x < —y, we have V(z,7T) = 0. 
Theorem 3.3 implies that the truncation error of the shifted approximation is O(h). 


This implies that the implicit method is consistent. 


Using the above approximation, we have that for 2 = 1,---,m—1andn = 
eee oe 
yor _yn yr = yo 1 a 
a a a i- | 


Sy ey) 


nt+1 
" ha Sa 
k 


=0 


AT h 


After re-expressing the above equality, this creates a system of linear equations, 


which can be rewritten as the following matrix equation: 
Av™tt =v" + Atk? 
where V4 = [pyre 
V" + ArF” = (0, (1 — yAr)V;", (1 — yAT)VZ, -.., V(y, 7)" 


and A = [A;,;] is the coefficient matrix, given as follows: 


0 j>it2 

082 pes 
Aga) 1 St- 952 Gat 

— 57 — got i ee a 

~ 94-541 9S je | 
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Let » be the eigenvalue of matrix A, so that for some nonzero vector X, we have 


AX = \X. Choose index i so that |x;| is maximized where x; represents the ith entry 


in X for 2 =0,...,m. Therefore, we have S- Aj,j;2; = Ax;, which implies, 
j=0 


j=0;5A% 4 
Thus, 
ce. 1+ 42 (1 = “1 —_“ ln b oye ora | i# {0,m} 
1 i= {0,m} 
LG 
Since |—4| < 1 and g; > 0 for any j = 0,2,..., we have 
Xi 
i+] sh itl 
Ss Gi-j+1 —l< S- Gi-jHi S 1 
j=0,54i Vil 5-05 Hi 
and, 
i+1 i 
mt > gi-j41|| <0 
j=0;jA1 Ui 


This implies that all eigenvalues of A satisfy |A] > 1. A is an invertible matrix so 
A™! exists. The eigenvalues of A~! satisfy || < 1 and therefore, the spectral radius 
of A~! satisfies p(A~!) < 1. Note that the error vector €° of V®° results in the error 
vector ¢! of V', related by et = A~‘'e®. This implies that |Je*]| < ||e°||. The method 


is unconditionally stable and the theorem is proven. 


In the next section, we will apply the Finite Difference Method using the shifted 
Grunwald approximation of the fractional derivative. We will mimic the mathe- 
matical contents presented in Meerschaert and Tadjeran for European vanilla option 


pricing. [9] 
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3.3. STANDARD FINITE DIFFERENCE METHOD 


2 
Let h = ala x; = —p+ith for 1 = 0,---,m, let t, = nAr where Az represents the 
m 
1 
incremental change of time in the interval 0 <7 < aca sec ST for n = 0,--- ,m, 
and: Vi = V(gecr,,). 


The shifted Grunwald approximation for the fractional derivative is given by the 


formula: 
O°V (az, T) 1 — 1 ST(k-a) 

= l V(a —(k—1)h,rT). 

axe T(—a) Rae he > T(k+1) ae Jy) 
Agai ill let Beso) be the normalized Grunwald weights. Let 

in, we wi = 
a Aras w= (a) ha - 

s(z,T) := (y —1)—— —)V be the source/sink term of the Black-Scholes Equation. 


Ox 
Under the interval —y < x < y, the finite difference method implies that 


Vi, -V;" 
sf = (ti, Ta) = (Y — EL — P+ O(h). 


Therefore, by using the shifted Grunwald approximation, we see that the Black- 
Scholes Equation under the finite difference method is given as: 


yor _ Vr 1 i+1 
a t VA ob 3”, 
Ar he dk i—k+1 a 


Solving for V,"*" yields that for —y < x < y, we have 


Ar At Fa a 
Vert = —goVii + (1+ —a] V+ — YS 9eV ee eg + sp Ar 
he he A ae 
2 
with h = ate and s defined above. For 7 = 1,---,m-—1 and for each n = 
m 
1,---,m—1, the m — 1 equalities degenerate into a matrix equation given below: 
Va I Jo 0 0 Ve 92 
as ge on Jo 0 vy 93 
AT Ar n > 
Yee ae pa | 98 g2 Ge. Ge 7) V3 | + Fa V0 ga| + ATSn 
yo Gm-1 Gm-2 Ym-3 J Vy Im 


with 


1-y-hy y—1 0 0 Vi 0 
0 1-—y-hy y-1 0 Vy 0 
os» ch y—-1 
SS > 0 0 1—y—hy 0 |= 0 
h h 
0 0 0 + Ll-y-hy Aba 1 


Solving the above matrix equation using Gaussian Elimination requires a compu- 
tational cost of O(m?) operations and a storage cost of O(m?) per time step. Note 


that for all n = 0,--- ,m, 


oe 0 for European call option 
(0) —S 
Keo for European put option 
and, 
i er for European call option 


0 for European put option 


Also, note that for all 7 = 0,--- ,m, 


yo max{e™ — k,0} for European call option 


max{Kk — e*, 0} for European put option 


In the next chapter, we will present an efficient finite difference method to solve 
the Black-Scholes equation under FMLS model. We will attempt to reduce the stor- 
age and computational cost of the method. In doing so, we must add a necessary 
restriction to the initial conditions inconveniently not satisfied by the problem as 


stated. 
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CHAPTER 4 


FAST NUMERICAL METHOD 


Under the standard finite difference method, recall that the price of a European 


vanilla option V(z,7) satisfies the following recursion: for —p < x < y, we have 


AT AT TAG esta 
ver = ho 901 as ( fr he ws) ed te ho ys OV" as sp Ar 
k=2 
with h = ae and, 
m 


The above recursion creates a system of equations which can be re-expressed as a 
matrix equation. Under the standard finite difference method, the matrix equation 
requires a storage of O(m?) and a computational cost of O(m*). In this section, we 
propose a more efficient method by reducing the computational cost in the inversion 
of the coefficient matrix. The result of this inversion method reduces the storage 
cost from O(m?) to O(m) and reduces the computational cost to O(mlog m) at each 
iteration. 

In this chapter, we will discuss about the limitations of the direct O(m log? m) 
finite difference method and propose a fast locally conservative finite volume method 
that similarly reduces the storage cost and computational cost of the standard finite 
difference method. In using this method, we require to re-express our equation in a 


conservative way. 


59 


4.1 LIMITATIONS OF UsING DirRECT O(m log? m) FINITE DIFFERENCE METHOD 


Let V(a,t) be the price of the European vanilla option with respect to log price x and 
time t. Let r be the risk-free rate of interest, ¢ > 0 be the volatility, and 1 <a<2be 
the fractional parameter. Recall that after a change of variables, for —oo < x < o, 


the FMLS Fractional Black-Scholes Equation is given as follows: 


with boundary conditions V(x, 0) = II(x) for H(z) representing the payoff function 
ar 


for European vanilla options defined previously and y = -2 COs >. 

For a direct O(m log? m) finite difference method approach, we require the extreme 
values of the log price to take equal values. [13] Note that for maturity 7 and log price 
x, the price of a European call option V(x,7) under the FMLS model is monotonic. 
Specifically, recall the analytic solution for European call options under the FMLS 
model: for log price of a European call option x and time T = 50" sec ats —t) 
satisfying 0 < t < 7, the price of the European call option is analytically given as 


follows: 


dy 


dy 
V,(a,7) = e | exp (—1 —1/*k) fag (|hl) dk — ket fa.o (|kl) dk 


—log kK —(1— 2 
where ie oe 7 Cie and -y = —-- cos — and 
FALo oe 2 
1 ( = +,+) (3,4) 
fo.o (\k|) = Hx ||h| i 
(0, 1) (x5) 


when ASS is the Fox H-function. 
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Taking the derivative with respect to x yields: 


OV. x 1/a 1 x at 1/a 
ae = ev exp (—1 -— 7 /d:) faalldal zs +€ ie exp (—r — r/R) faso (Ikl) dk 
1 
— Ke” fa olla) a 
- 
dy 
= | exp (-r — rk) foo (\k|) dk > 0. 


By Theorem 2.4, we have, 


dx Ox 
Or, 
OV, 2 [~ 1/a 
eo te ip exp (-r-7 k) fo (\k|) dk < 0. 


Therefore, we have shown that if V represents the price of a European call op- 
tion, V(x,7) is monotonically increasing with respect to log price x. Similarly, if V 
represents the price of a European put option, V(2,7) is monotonically decreasing 
with respect to log price x. 

However, it is not possible to truncate the tails so that there exist two endpoints 
with equal values. If there exists LZ and R so that V(L,7) = V(R,7T), by Rolle’s 
Theorem, there exists a relative minimum or maximum in the interval L < x < R. 
This is a contradiction since the European option value function is a monotonic 
function. 

A naive approach to circumvent this issue is to study the Taylor series represen- 
tation of V(z,7). This will require the Taylor series representation for Levy stable 
density functions. In general, the Levy stable distribution can be expressed as the 


real part of a simpler integral given below: 


F(eie,8,6,n) = bRe| [exp (it(e — n) — (e* (1 -i9tan %))] 


By expressing exp (-(ct) (1 — 76 tan on)) as a Taylor series and reversing the 


order of integration and summation yields, 
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I, 
f(x0; 86,4) = -- 
a xr — [lb 


ap oe" : ) Tan + 1) 


which is valid for all 2 # where q := c*(1 —i8 tan $). Thus, 


1 i ant+1 
fao(@) = f(z; a,0,1,0) = —Re se mc ) T(an+ 7) 
a n=1 = zt 
OF, 
1< (=1)"F (an 1) —an—1 antl 
foo(X) = a 2 al R i 
Note that, 
R Cage ne ( uf + isin ae = cos (= + ) = —sin ao 
e |2 = Re || cos 5 + tsin 5 = 5 5) = a 
Thus, 
1S (-1)""T(an+1) . (/man\ _.,-1 
Ffao(X) a S- al sin ( 9 )a ‘ 


For general Levy stable density function, we have the following asymptotic: 


Fic, 8, 6,4) ~ —e(1 + sen(2)8) sin E(a + fel. 


Therefore, 


1 
fao(x) ~ —sin —I(a+1)|2|7-'°. 
1 


However, in this naive approach, there are many limitations in the direct O(m log? m) 
method. In performing numerical experiments using the Taylor series representation 
of the closed-form solution V,(x,7), we cannot express the integrals in terms of stan- 
dard built-in functions. Due to the series conditional convergence, the question of 
the accuracy of the numerical approximation also arises. 

The low convergence rate at « — oo of the Levy density function f(x) has 


the added difficulty of computing semi-infinite integrals in our closed form using any 
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numerical methods. Although the generalized Laguerre-Gauss quadrature efficiently 
calculates integrals of this form, a naive truncation proposed in the Taylor series ex- 
pansion, specifically the infinite series of the Grunwald-Letnikov fractional derivative, 
yields a finite difference method that is unconditionally unstable. The infinite series 
also does not accurately approximate put values for all x satisfying « < Ke%-77 
since the assumption of d, > 0 is necessary to define the Taylor series expansion of 
the Levy density function fo o(x) for all 1 <a < 2. 

To avoid these issues, we will alternatively use fast finite volume methods which 
similarly use banded coefficient matrices to reduce the computational cost from O(m?) 
to O(mlogm) and storage cost from O(m?) to O(m). For the application of this 
method, we relieve the restrictive boundary conditions of u(L,7) = u(R,r) = 0 and 
re-express the Black-Scholes equation in a conservative way. This calls for a crucial 
definition and some additional discussion. More information can be found in Cheng, 


Wang, and Wang. [6] 


4.2 PRELIMINARIES FOR FAST FINITE VOLUME METHODS 


Let a= 2-8 for0< 6 <1. The left and right Caputo fractional derivative of 


order 1 — £ are respectively defined as follows: 


DADg2) = Fr f (e— 9)" *al(shds 
Dy Dale) = peg [ (s— 2) 9's) 


Recall that after a change of variables, for —oo < x < oo, the FMLS Fractional 


Black-Scholes Equation is given as follows: 


OV 


OV OV(a,T) 


| — VV 
Ox® Y 


with boundary conditions V(x, 0) = II(x) for H(z) representing the payoff function 
ous 


2 
for European vanilla options defined previously and y = — =" cos SF 
o® 
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Let V(x,7) = exp(—77T)u(z,7). Then, we have the following: 


oe = exp(—yr) Se — yuexp(—17) 
Bp = exP(—IT) a- — yuexp(—IT 


Furthermore, we have, 


O°V ce T(a + 1)D2-"u(a, 7) D? exp(—yr) 
One eC, a T(a—n+1)n! 


= Dtu(xz,7T) exp(— yr). 


Details of this proof can be referred to Osler. [10] Our partial differential equation 


becomes the following: 


Or equivalently, 


_ Ou 
ae pp De Putt) = =4) 5 


The above variant of the Black-Scholes equation under the FMLS model is re- 
expressed in a form renowned as divergence form. ‘This is the required form 
necessary to successfully apply the fast finite volume method. Unlike the direct 
O(mlog? m) finite difference method, we do not require u(—y,T) = u(y,T) = 0. 

Continuing our construction for this method, let us define a temporal partition 
on the interval 0 <7 < —50° sec — bye = Aros A... wih ar = 
—0° sec > We will use the finite difference approximation of the time derivative 


to re-express the partial differential equation at time 7, as follows: 


u(Z,T) — U(£,Tm-1) 9 
At Ox 


(.D5°u' (a, ™n)) =f (a) 


Ou 
where f(2,7) := (y — 1)—(2, T) represents the source/sink term. Equivalently, 


Ox 
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U(Z, Tm) — Aro (Gerace. qn) S ATT (tte) Ue, Tea) 


Next, define a uniform spatial partition on the truncated interval of the log price 


2 
[—-y, v] by letting 2; = th for i = 0,1,....m with h = ae Furthermore, define 
m 


L,-1 = (aj-1 + 2;)/2 for i = 1,...,m to be the midpoint of the interval [z,_1, z;]. 


9 
Then, by integrating the governing equation on (2,1, 441) for i = 1,...,m yields the 


following: 


/ a u(x, ™m)da—Ar ,D5?u' (x, Tr) 


x. ved 4 reed 0 
oS Ar f eee f(a, tm)dx-+ f "2 ula, Tr-1)dx. 
ae Cad asl 


Let S;,(—y, y) denote the space of continuous and piecewise-linear functions with 
respect to the spatial partition that vanishes at x = —y and x = y. Define the nodal 


basis function ¢;(x) for k =1,...,m— 1 as follows: 


Lea, —2£ 
be(x) = as x E (Le, Le41| 
0 elsewhere 


Let uZ := u(x, T). Then, ua(x, Tm) € Sr(—y, y) can be represented as: 
m-1 
i (25-7) = — UpPr(z). 


k=1 


Thus, for all i = 1,...,m—1, the finite volume scheme is formulated as: 


oe m—-1 Eid 
= Ar | ? f(x, T)dx + De ay *? 6,(a)dx. 
oe j=l 8 


Let ff := f(x,,T) and define the vectors u® = [uf,u%,...,u™_,]? and f" = 
fr, fe, f"% |" be the vectors corresponding to the numerical approximation and 


the source/sink terms at time step 7,. Let M = [Mj,]/"—) be the mass matrix and 
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B” = [BP 


Ha 


;;-1 be the stiffness matrix at time step 7,. The entries of the matrices 


M and B” and vector f” are given as follows; for 1 <7i,7 <m-—1: 


B= i. “2 D,.D5° Dé, («)dx 
2 


ie CR 


1 
2 
We will shortly determine the entries for and B”. The entries of vector f” are 
simplified as follows: 
i Tih _ os 74d _ = non 
f= f"* Fe ta)de = (71) unstable = (7-1) (why hg): 


The finite volume scheme can be re-expressed in the following matrix equation: 
(M — ArB")u" = Mu” +Arf” 


The mass matrix M takes the same form as its counterpart for the classical dif- 


fusion equation. It is given as follows: 


le 1 0 0 ol 
16 1 0 
1/0 1 6 
M=-h 
0 
Oa, Faye HES 
OP 0 Sek FOr “76 


The difference between the finite volume scheme and the classical diffusion equa- 
tion lies in the properties of stiffness matrix B”. By observing the support of the 
differential operators, we can conclude that the stiffness matrix is a full matrix which 
requires a storage cost of O(m?) and a computational cost of O(m?). Many meth- 
ods, including the finite volume method being discussed in this chapter, have been 


proposed to effectively reduce these costs. 
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Under algebraic manipulations, one can deduce that the entries of the stiffness 


matrix can be rewritten as follows: 


l (3) 
Pies ane a) 
I T(B + Dai-BF-3 


where g? ) ig defined as follows: 


Jj 


Alternatively, we can decompose the stiffness matrix as follows: 


1 
"= Pep Dare? 
where, 
9 GF Gs Gon 
Ga Ga Sat ae ees 
(ee 
dn om GP of 
Gee dees ee a 2g 


Note that B” does not depend on n. We will rename B := B” to emphasize this 
relationship between the stiffness matrix and the time step n. Our matrix equation 
now becomes: 

1 1 
—_M—~B a —_M n-1 n 
(a ) e Ar oa 


Note that we can conveniently compute aM as follows: 
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1 6 1 0 
AT ~ 962-87 ~* 2 0 
0 6 1 
0 O -:-: O 1 6 


Contrary to the standard requirement of O(m?) storage cost for any full matrix, 


the matrix B requires O(m) storage. 


4.3 Fast O(mlogm) ALGORITHM FOR THE EVALUATION OF Bu 


In this section, we discuss a more efficient method of evaluating the matrix-vector 
multiplication of the stiffness matrix B and some vector u. Firstly, define the following 


matrix G as follows: 


O. Gar ge 
Ge (0 a. 08 
= 
ae Os 293 
Ge ig aarp: . 10 


and define the following (2m — 2) x (2m — 2) circulant matrix as follows: 
G 

Com—2 an ee 

G 

The circulant matrix C,,-2 can be decomposed as follows 


Cia Fyy,—odiag (Fam —20) Fom—2 


where c represents the first column vector of Cgn~2 and Fyn—2 represents the 


(2m — 2) x (2m — 2) discrete Fourier transform matrix with entries given by 
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1 ie mgd 
V 24 F em 


Given this, we describe the procedures necessary to evaluate this matrix-vector 


Fym_2(9, 1) = 0<7,1<2m—-3 


multiplication efficiently. 
1. Rewrite u as Uzm—2 = [u?,0]” so that, 


Gu 


Com—2V2m—2 = 


Gu 


2. Similarly, evaluate the matrix-vector product Wan_2 = F2m—2U2m—2 in O(m log m) 
operations by using fast Fourier transform. w2,_2 is the discrete Fourier trans- 


form of Uam_2. 


3. Now, we evaluate the matrix-vector product Vam_—2 = Fam—2C2m—2 in O(m log m) 


operations by using fast Fourier transform. 


4. Next, evaluate the Hadamard product zam—2 = Wom—2'Vam—2 = [W1V1,°°* , Wan—2Van—2|" 


in O(N) operations. 


5. Lastly, evaluate Yyom—2 = Fop_9%2m—2 in O(mlogm) operations using inverse 
fast Fourier transform. This yields: 
Gu 


Yom—2 = Com—2V2m—2 = 


Gu 


6. Multiplying by the scalar to the vector y2m—2 in O(m) operations 


1 
T(8 +1)h'-8 
yields Bu. 
By the above procedures, it can be deduced that the computational cost of evalu- 
ating Bu for stiffness matrix B and some vector u is O(mlogm). This is a significant 
improvement compared to the standard finite difference method’s computational cost 


of O(m) operations. 
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4.4 FAST CONJUGATE GRADIENT SQUARE METHOD WITH O(m) STORAGE 


Since the stiffness matrix B is full, the direct method of computing the matrix- 
vector multiplication requires O(m?) per time step. Given that we have deduced an 
algorithm to efficiently calculate Bu, it motivates us to consider the conjugate gradi- 
ent square method (CGS), a generalized method of the standard conjugate gradient 
method. This will allow us to accelerate the performance of the finite volume scheme. 

We avoid using the standard conjugate gradient method as this method does not 
apply for nonsymmetric systems The residual vectors cannot be made orthogonal with 
short recurrences. The pseudocode for CGS on the finite volume scheme discussed 


throughout this chapter is provided below: 


Algorithm 1 Conjugate Gradient Squared Method 


1: procedure AT EACH TIME STEP T,, CHOOSE uw =u"! AND COMPUTE 1 = 
f” — A"u. CHOOSE # (FOR EXAMPLE, # = 7) 


23 for 1=1,2,--- do 

3: Pi-1= FFG) 

4: if p;_; = 0, the method fails 

5: if 7 = 1 then 

6: wh) = 

ui pos yo 

8: else 

9: Bi-1 = pi-1/pi-2 

10: w= HN +B g&—) 

11: p® = w + Bi 1(g + Bip) 
12: — A”yl 

13: Qi = Pi- 6 

14: gq? = w® — ag 

15: u =ulY + aj(w + q) 

16: g = A”(w + q®) 

17: re) = 7 ace 

18: 5 = ||f" — Aru|| 

19: check convergence; continue if necessary 
20: u® =u) 


The above pseudocode requires O(m?) operations to compute the matrix-vector 


multiplication. By the previous section, it has been shown that the computational 
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cost can be reduced to O(m log m) operations. Applying the fast algorithm discussed 
in Section 4.3 onto FCGS method yields a fast finite volume method that efficiently 


computes the numerical values of the price of European vanilla options. 
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CHAPTER 5 


NUMERICAL EXPERIMENTS AND OBSERVATIONS 


In this chapter, we will carry out numerical experiments to investigate the perfor- 
mance of Euler’s backward method for the finite volume scheme discussed in Chapter 
4. We will also implement the Fast Conjugate Gradient Squared Method to acceler- 
ate the performance of the finite volume scheme and to observe its reduction of the 
CPU time. Lastly, we will simulate the Black-Scholes model under the FMLS process 
for European puts with different choices of a to emphasize the effects of the return 


distribution when compared to the underlying price. 


5.1 SIMULATION OF S&P 500 OPTIONS MARKET 


Firstly, consider the one year time interval between September 15, 2014 to September 
15, 2015. Using recently observed data from the S&P 500 options market in this time 
interval, we note that Chicago Board Options Exchange Volatility Index (VIX), a 
standard measure of implied volatility, consistently fluctuates between o = 0.1088 and 
ao = 0.4074. For the purposes of our numerical experiments, we will let ¢ = 0.1601822 
represent the average daily VIX observed in this time frame. The data can be observed 
and retrieved from YAHOO! Finance. 

Since the 2007 financial crisis, the U.S. Federal Reserve have kept interest rates 
between 0% and 0.25% in attempts to stabilize the U.S. economy and the financial 
system. This consequentially creates unnecessary inflation on many European vanilla 
options. Recent trends have suggest the U.S. Federal Reserve to consider raising 


interest rates to reduce inflation. For the sake of our numerical experiments, we will 
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BS Pricing Formula for European Puts Under FMLS Model 


Option price 


5 6 7 8 9 10 11 12 13 14 15 
Underlying Price 


Figure 5.1: European puts at different a@ values. Model parameters are kK = $10, 
r= 0,025, and 7 = 1, 


let r = 0.0025. 

Lastly, we will assume that the maximum asset price is $150. This suggests that 
we will truncate our log price interval to the following interval for both finite difference 
scheme and the finite volume scheme: —5 < x < 5. For our numerical experiments, 
we let p = 5. 

Under these realistic conditions, we are able to simulate the Black-Scholes equation 
under the FMLS model for varying parameters a. Provided above depicts the overall 
effect of European puts under different choices of a assuming that strike price kK = 
$10. In our simulation, we choose a = {1.25, 1.5, 1.75, 1.875,1.95}. As expected, 
as a — 2, the tail behavior of the return distribution becomes flat and the effect 
of leptokurtosis diminishes. Under put-call parity, we should also expect a similar 
behavior for prices of European call options. 

We will also discuss the CPU time usage for the methods presented in this thesis. 


Under the finite volume scheme (FV), we consider a spatial and temporal partition 
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defined by the midpoints of consecutive nodes and solve the Black-Scholes equation 
under FMLS under these points using Gaussian Elimination. Furthermore, we will 
implement a Fast Conjugate Gradient Squared acceleration (FCGS) onto the finite 
volume scheme to reduce the CPU time. The results are compiled in the following 


table. 


Table 5.1: CPU Time Usage of Numerical Methods. 


m | FV with GE | FV with FCGS 
23 0.094 s 0.016 s 
oe 0.4414 s 0.153 s 
on 1.706 s 0.478 s 
sid T2388 1.235 s 


As noted in Table 5.1, it is expected that the finite volume scheme with Gaussian 
elimination requires more CPU time to successfully determine the prices of European 
vanilla options under any choice of a. Using various choices of m, we can observe 
that the times increase significantly under this method. Recall that we require a 
computational cost of O(m?3) and a storage cost of O(m?) operations. This suggests 
that the CPU time increases exponentially as we scale m as observed in the table 
above. 

When applying the Fast Conjugate Gradient Squared method, we successfully 
reduced the computational cost to respectively O(mlogm) operations per iterate. 
Thus, we should expect that the CPU time reduces for all choices of m. Again, this 
can be observed by the results compiled in the table above. Under this acceleration, 
we reduce the storage cost from O(m?) to O(m). We should expect a reduction of 
CPU time since the finite volume scheme under this method does not require much 


storage. This again can be observed in the data presented in Table 5.1. 
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We now turn to the accuracy of the numerical methods. We will compute the 


relative error in the Lz and L., norm. Recall that Lz error is computed as: 
1/2 
6 = (= lV = ve) 
And the ZL. error is computed as: 
Eoo = max |V — V;| 


where V represents the actual price of the European vanilla option under fixed 
qa and V; represents the numerical value evaluated at x; for T = 1. The results are 


presented below: 


Table 5.2: Accuracy of Numerical Methods 


m €9 for FV €oo for FV 


2? | £0531 %10-* | 3.1012 x10-3 


210" bAbS2 X 105 | 1.596% 10 


2° || D890 e102" | 71088 10. * 


22) "OF e103 | 3.2579 0 


As expected, we note that as m increases, €2, €,, — 0, since we require finer tempo- 
ral and spatial partitions in the simulation of the Black-Scholes pricing formula under 
the FMLS process. Also, note that under the finite volume scheme, our simulation 
may motivate more methods to efficiently compute the desired option prices given 


parameter a since these errors are insignificant to the application of option pricing. 


5.2. CONCLUDING REMARKS 


In this thesis, we have successfully analyzed the behavior of the Black-Scholes model 


under the FMLS process for varied choices of a, determined a closed-form solution 
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to represent a variant of the Black-Scholes pricing formula in terms of a, and used 
numerical methods to simulate the recent behavior of the S&P 500 options market ina 
one-year time frame. We have also considered finite volume schemes to further reduce 
the CPU time by decreasing the computational and storage costs of the standard 
finite difference method, using banded coefficient matrices to efficiently compute the 
solution of the system. Lastly, we performed an applicable simulation emphasizing 
the efficiency of the finite volume method and have successfully achieved our desired 


results. 


2 
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APPENDIX A 


DIFFUSION MODELS AND RANDOM WALKS 


In this Appendix, we provide a mathematically justified proof of the diffusion equation 
using applications in random walks. We will begin by proposing the Central Limit 


Theorem. We will provide no proof of the Central Limit Theorem in this Appendix. 


Theorem A.1. (Central Limit Theorem) Let X1,--- , Xn be independent, identically 


distributed random variables satisfying E(X;) = m and Var(X;) = 0? > 0. Let 


Sy t= Xi +---+Xy. Then fora,be R,a< b, 


=s b 
im, P (as A x) =f’ Fav 
noo Jno In Ja 


In our mathematical justification of the diffusion model, we require the following 


theorem. We will provide a proof of the theorem. 


Theorem A.2. (Laplace DeMoivre Theorem) Let X1,--- ,Xn be independent, iden- 


tically distributed random variables satisfying, 


P(X;=1) =p and P(X; = 0) =4, where i € {1,--- ,n} 


for p,q=>0 andp+q=1. Define S, := X,+---Xn. Then fora,b€ R, a < b, 


Sy — np 1 bg? 
lim P{a< re eee i > de. 
na¥00 («< /Npq J2Qr a ‘ is 
Remark: One can easily compute that E(S,,) = np and Var(S,,) = npg. The 


Laplace-DeMoivre Theorem is a specific case of the Central Limit Theorem. 


75 


Sn. — Tp 


—n 
Proof. Define S) := to be a random variable with value x, = E 


Vrpq V pd 
probability p,(k) = (() p*q”*. Then, for a < b, 


and 


n n) 
Sgt = kan—k 
Passi <b)= D0 pnalk > (i): q ye ie =a pe q’ 


ac<ap<b a<ap<b a<ap<b ( 


By applying Stirling’s Formula, which states that as n — oo: 
n! =e "n"V/2mn(1 + o(1)) 
we have that as n > oH, 


en” /2nnp*qr-* 
Plax Si <b)= (1 + o(1)) 
oy ekkky/2nke~-*) (n — k)"-*,/2n(n — k) 


7 oe : rE k) (2) ee (1 + o(1)). 


Now we propose the following lemma: 


Lemma A.3. We have the following: 


R= 
ae Note that, 


npd 
Iq Iq (2) k; 
np np \ ./npq np 
it: | Dp = | ‘p zi vn) n—k 
nq 


Given that as y > 0, we have log(1+ y) = +y—- & ae O(y’). Therefore, we have 


Proof. Let x = x, = 


and, 


—(np + 2./mpq) (ee 2 + O(n-?) 
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And similarly, 


in oe nk 


= —(nq — x./npq) (- Renae “) + O(n-2) 


nq 2nq 


Adding the expressions and simplifying yields, 


k n-k 2 
i, tel) (yo 
N03 Fang Pe n— 2 


7% 


This implies the statement of the lemma. The lemma is proven. 


k —np 


Vp 


Again, let x = 7, = . Therefore k = np+2,/npq and n—k = nq—x.,/npq. 


Thus, 


n 1 
k(n—k) npg 


Continuing from our calculuations, as n — oo, we have: 


(1+ o(1)). 


1 1 a 
PGs 5) <5) =—— ——e 7 (1 1)). 
(a =e Drage Oto) 


Note that the right hand side is simply a Riemann sum approximation as n — oo 


of the following integral: 


This completes the proof of the theorem. 


We will now proceed to derive the traditional diffusion model using random walks. 
Given spacing Ax > 0 and time duration At > 0, define the following two-dimensional 


lattice: 


{(mAz, nAt) |meZ,n Ee Nv {OF}. 


ws 


Consider a particle in starting position x = 0 at t = 0. We allow the particle to 
only move to the left and to the right. Let the probability of the particle moving left 


1 3a3 3 3 : 
5- Therefore, the probability of the particle moving right Ar 


an amount of Az be 
units is 5. 
Let p(m,n) be the probability of a particle at position mAz and nAt. Then, we 


have, 


Furthermore, we have 
1 1 
p(mn,n +1) = sp(m — 1,n) + 5pm + Ln) 
and equivalently, 


p(m,n +1) ~ plm.n) = 5 [p(m — 1,n) ~ 2p(m,n) + p(n — 1.) 


Ar)? 
Assume that ( ot) =o? > 0. Then we have, 


p(m,n+1)—pl(m,n) — 0? p(m—1,n) — 2p(m,n) + p(m — 1,n) 


At 2 (Ax)? 


Letting At > 0, Av > 0, mAz > z, and nAt - t, we see that p(m,n) > u(z,t) 
where u(x,t) is the probability density function of a particle at position x and time 
t. We also see that by passing limits, the difference equation becomes the traditional 


diffusion model, i.e. 


Ou(x,t) 0? u(x,t) 
Ot 2 Orr ° 


Now we will solve the traditional diffusion model by using another interpretation 
of random walks. Again, consider the particle with the same conditions previously 


mentioned. Let X(t) be the position of the particle at time t = nAt, and define 
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le 1 
——— S° Xi where X; are independent random variables satisfying P(X; = 0) 5 
i=1 
ih 
Z 
S, is a random walk that determines the number of moves to the right at time 


t = nAt. Therefore we have, 


1 
and P(Xe= 1) = 5 Note that Var(X;) = 


X(t) = S,Az + (n — S,)(—Az) = (28, — n)Az. 


Note that, 


Var(X(t)) = (Az)*Var[(2S,, — n)] 


= 4(Ar)*Var[S,] = 4n(Az)?Var[X|] 


Continuing from our calculations, we have 


g 8 g 2 
X() = 54 mae ( "2 | /nAg = | “2 | Vort. 
ee <2 


Applying Theorem A.2 yields the following result: 


This implies that the probability density function u(x,t) satisfies 


i} ig 
tet) = are 353 | 


The above solution also solves the traditional diffusion model. 
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APPENDIX B 


STOCHASTIC CALCULUS 


We will begin Appendix B with a few necessary definitions that play a crucial role in 


the derivation of the Black-Scholes model. 


Definition B.1. A real-valued stochastic process W(-) defined on some probability 
space (Q,U, P), is called Brownian motion (or Wiener process) if W(-) satisfies 


the following conditions: 


e W(0) =0 almost surely, 
e W(t) —-W(s) is N(0,t—s) for allt >s>0, 
e for all times 0 < ty) < tg < ... < tn, the random variables W(t,),W (te) — 


W (t1),...,;W(tn) — W(tn_1) are independent. 


The stochastic process W is imperative in reflecting how stock option market 
pricing behaves under the most simple conditions; namely, it attributes to the behav- 
ior’s resemblance of a random walk which inspires the construction of the traditional 


diffusion model. We will derive the diffusion model later in this section. 


Definition B.2. Let W be a Wiener process. For times0 <t <T’, we define Ito’s 


Integral to be the following: 


WT) TT 
-5. 


T 
i Wane 
0 2 


We will assume that all stochastic integrals, including the Ito’s Integral, exists and 
are well defined throughout this paper. More details on the existence of stochastic 


integrals can be referred to Evans. [?Evans13] We propose the following definition: 
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Definition B.3. Let X(-) be a real-valued stochastic process satisfying 
X(r) = X(s)+ f pat + f odW 
where 4s and o are real-valued progressively measurable processes satisfying 


, - 
D (/ wit) <0o and (/ oat) <0 
0 0 


and s andr are times satisfyingO<s<r<T. Then for0<t<T, we say that 


1c] 


X(-) an Ito process with stochastic differential dX = dt + odW for drift yu 


and volatility o. 
We will begin this section by proposing Ito’s Chain Rule. 


Theorem B.1. (Ito’s Chain Rule) Let f(x,t) be a smooth function and let X(t) be 
an Ito process with stochastic differential dX = dt +o0dW. Then Y(t) := f(X(t),t) 


is also an Ito process with stochastic differential: 


Or equivalently, 


2 
dY = (F gee i pe f) dt 4 gla. 
Ot bi 


Ito’s Chain Rule plays a crucial role in the stochastic derivation of the Black- 
Scholes Equation. In order to prove Theorem B.1, we need to propose and prove Ito’s 


Product Rule. 


Theorem B.2. (Ito’s Product Rule) Let X(t) and X2(t) be Ito processes satisfying 
dX, = pdt +o0,dW and dX2 = podt + oodW. Suppose that p14, 2,01, and og are 


real-valued progressively measureable processes satisfying, 
a T a4 iP 

D i \uildt| , E | \tua|dt | < co andE | ofdt|) , E i: a3dt | < oo 
0 0 0 0 
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for0<t<T. Then, X,X2(t) is an Ito process satisfying 
d(X,X2) => XodX1 + X1dXo + 1 02dt. 


Remark: Note that the integrated version of Ito’s product rule yields Ito integration- 


by-parts formula: 
| " XodX = Xy(r)Xo(r) — X1(s)Xo(s) — i X1dX> — , Gide 


Proof. Choose 0 < r < T. Let F(t) := U(W(s), Xo) be the o—algebra generated 
by Xo and W(s) for 0 < s < t, where Xo is a random variable independent of the 
future of Brownian motion beyond time t = 0. Firstly, assume for simplicity that 
X1(0) = X2(0) = 0, pi(t) = pa, peo(t) = pe, o1(t) = o1, and oo(t) = oo, where 
/11, [42,01, and o» are time independent, *(0)—measurable random variables. Then 
for t > 0, we have X,(t) = pyt + 0,W(t) and X(t) = pt + ooW(t). 

Thus, 


[ OR AR Ride eoxeadt = if ” Ridin + Noite i: “Spt Roo aw [ “erat 
= [at + 0,W) po + (pot + 0oW) pdt 
+ [at + 0,W)o2 + (Mat + o2W)o,dW + ip o\02dt 
= py per? + (orplg + o2}11) ([ Wdt + if idW ) 
+ 20109 a WdW + ojoor. 
By Definition B.2, we have af WdW = W?(r) —r. We propose and prove the 
following lemma: 
Lemma B.3. For r > 0, a tdW + _ Wadt =rWi(r). 


Proof. Let P” = {0 = tj < t] <--- <i, =r} be a sequence of partitions of the 
interval [0,7] with |P"”| + 0 as n — oo. Therefore, by definition of the stochastic 


integral using Riemann sum approximation, we have: 
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Mn -l 


[ taw = kim, > RCW (EE) — WED). 
k=0 


Similarly, since t+ W(t) is continuous almost surely, we have, 


r Mn -l 
[Wet = jim OS Wt) - 0. 
k=0 
Therefore, we have, 


Mn-1 


: tdW + is Wat = lim S° tf(W (tear) — W(ER)) + W (thy) (thar — #2) 
k=0 


Mn—-1 


= jim ~ thy1W (thi) — teW (te) 
k=0 


WV (tn) — toW (to) 


With these identities, we have 


i. Xod X41 + X1dXo + o109dt = [i fler? + (O1pl2 + Oop, )rW (r) + o102W?(r) 


= X1(r)Xo(r) 


which is equivalent to the Ito integration-by-parts formula for the specific case 
of s = 0, X1(0) = X2(0) = 0, and py, W2,01 and og are time-independent random 
variables. A similar proof can be extended for the case of s > 0, Xi1(s) and X2(s) 
are arbitrary, and /11, 2,0, and 2 are constant F(s)-measurable random variables. 
Ito’s Product Rule holds for constant random variables. 

Suppose that [11, /42,01, and a2 are step processes. For each subinterval [t,, t,+1) 
for which 11, 42,01, and g2 are constant random variables, we can apply the previous 
steps and add the resulting integrals. Ito’s Product Rule holds for step processes. 

Now suppose that /11, [l2, 01, and a2 are general processes. We can select a sequence 


of step processes pi], 5, 07, and of such that 
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r ir 
e([ iat ~ maa) >and B { [ it~ pal] 0 
0 0 


and 


E ( [er = oat) — 0 and EF ( [es - oat] — 0 


0 


holds. Define 
t t 
X(t) = X,(0) + | yrds + i otdW 
0 0 
t t 
X(t) = X9(0) + | ynds + | ondW 
0 0) 


Note that Xj’ and X3’ have stochastic differentials dX?) = pdt + o/dW and 


dX} = phdt+oafdW. Since X/ and X3' are step processes, Ito’s product rule applies. 
This yields d(X? Xf) = X$dX? + X?dX} + ofotdt. Letting n — oo yields the 


statement of the theorem. 


Proof. Let X be an Ito process with stochastic differential dX = udt+adW. Firstly, 


consider the case f(x,t) = 2” form € N. We claim that 
1 
d(x™) = mX™ "dX + gi(m —1)xX™ 07 dt 


and prove the claim with induction on m. Letting m = 0,1 yields trivial cases 
and for m = 2, we have d(X?) = 2XdX +o7dt. This ie equivalent to the statement 


of Theorem 1.2 for X; = X» = X. Let us assume that 


1 
dx = (m _ TX + 5(m = 1)(m _ DX ede 


1 
= (m—1)X™*(udt + odW) + aim =1)(m—2)X"™ *o7 dt. 
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Then, using Theorem B.2, we have 


A, Gus teat) Gua 
= Xd(X™ 1) 4+ X™1dX + (m—1)X"™ *o7 dt 
ae ((m —1)X"aX + s(m TG 2)X"- Scat) 
+ X™ dX + (m—1)X™" 707 dt 


1 
=mX™ 1dX + gin(m —1)X™ 07 dt 


The induction is complete. We can generalize the above case for all polynomials 
f (x,t) in the variable x, because the stochastic differential operator is linear. Ito’s 
Chain Rule is therefore proven for all polynomials f(x,t) in the variable x. 


Now suppose that f(a,t) = fi(x)fo(t) for polynomials f; and fo. Then, 


d(f(X,t)) = d(fi(X) fa) 
= fi(X)dfo + fodfi(X) 


= fi(X) fidt + fr (F(x)ax i ; fll X)o?at) 


By the above calculations, Ito’s Chain Rule applies for polynomials of the form 
f(a,t) = fi(x)fo(t). This can be generalized for all polynomial functions f of the 
variables t and x since any f can be re-expressed as a linear combination of functions 
of the form f(x,t) = f(x) fo(t) and the stochastic differential operator is linear. 

Now let f(z,t) be a smooth function. There exists a sequence of polynomials f” 

aie Of of” oF ead Oye ony 


such that f” — f, ey a? Ox => Da ee) > Aa?’ By applying the 


previous steps, we can deduce that for all O<s <T, 


eOF n0y" Ps esl ane | s Ofr 
ra BW+ ort + | 3g Cow Bs: 


f"(s, X(s)) — f"(0, X(0)) = 


Letting n — co yields Ito’s Chain Rule for all smooth functions f(z, t). 
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APPENDIX C 


BLACK-SCHOLES FORMULA 


In this Appendix, we will transform the Black-Scholes equation into the traditional 
diffusion equation by making specific substitutions. We have determined the solution 
of the diffusion equation in Appendix A. Finally, we will use this solution to derive 
the Black-Scholes Formula mentioned in Chapter 1. 

Let C(S,t) be the price of a European call option with asset price S and time t. 


Recall that the Black-Scholes model satisfies the following differential equation: 


OC 0c 1 aC 


2 a2 - 
Lr H 50 8 a5 rc. 


and boundary conditions that for 0 <t < T and 0 < S < ow, we have C(S,t) ~ S 
as S + oo, and C(S,T) = max{S — K,0} for strike price K. 
2 
Firstly, let S = Ke*, t = T — —r, and C(S,t) = Kv(z,t). Then, we have 
a 


ge = log(S/K) and 7 = SOUT —t). 
O-. 07-0 ia 
OS OSdx SOx 


Pe aa 2 (52)- 12+12(2)- 10.18 


1 
Equivalently, we have dr = — 50 dt and . Thus, 


0S? dSOS OS\Sdr) ~—S?0x' SOS \dx) ~~ S? Ax" S? Ox?" 
Making the above substitutions to the Black-Scholes equation yields, 


2 
1 ,du | sae] + 50°"| 1 Ov | = 9aa| =" 


aa; ae ~~ §20n | §2 0a? 


Or equivalently, 


86 


2 
Redefining k := = and rearranging yields, 
a 


Under these transformations, the boundary condition of C(S,T) = max{S—k, 0} 
for strike price kK becomes v(x,0) = max{e* — 1,0}. 

We now make the following substitution to further transform the above differential 
equation into a diffusion equation. Let v(z,7) = yu(x,T) where y(x,7T) = exp(ax + 


Gr). Then, we have the following: 


Ov Ou 

OY = Bae, ryu + 9(0,7) 9" 

Ov Ou 

Ax = ay(x,T)u4 (eT) 5 

O?u ‘ Ou Ou 
aan avy(x,T)u + 2ay7(z, a, -- V2.7) 55 


Making this substitution to the above differential equation and dividing both sides 


by (z,7) yields, 


Equivalently, we have, 


Ou j Ou Pu 
7, = B+a talk—1)—k)ut (Qat+k—I)s + 5G. 


1 1 
We require that a = 5(l —k) and 6 = —4 +k)? to force the above partial 
differential equation to become a diffusion model. With these choices of parameters, 


we have, 


Ou Ou 


Or Ox? 
with wala) = ule. 0) = max {exp GG + 1)x) — exp iG - 1)x) Oh. 
The diffusion equation that we have derived in the paper and Appendix A all have 
boundary conditions independent of x. We generalize the solution of the diffusion 


equation by applying the following definition: 
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Definition C.1. Let u(x,t) satisfy the following diffusion model: 


Ou Ou 
Or Ox?” 
for some constant k and boundary condition u(x,0) = g(x). Then, Green’s 


function for the diffusion model is given as follows: 


iat) = 


— 4) 
_ dy. 
ac ze(-& akt ene 
Applying Definition C.1 for k = 1 and the change of variables y +> x + V2rTy, we 


see that 


Ue T ) 


=s= uo(a + V2Ty) exp (-© dy. 


Reverting the solution to the original variables yields the desired Black-Scholes 


formula given as follows: 


C(S,t) = SN(d\) — Ke"? N (dp) 


where, 
N(x) := ae exp (—5s?) ds 
And, 
1, _WealS/K) + r+ oT -1 
o/T —t 
Ane log(S/K) + (r — 50°)(T —t) _ PY eer: 


ovVvl—t 


More details can be referred to Dewynne, Howison, and Wilmott. [?Dewynne95] 
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APPENDIX D 


FOURIER TRANSFORMS OF DIFFUSION MODELS 


We will begin this Appendix by providing the statement of the Levy Continuity 


Theorem: 


Theorem D.1. (Levy Continuity Theorem) Let {X,}°2, be a sequence of random 
variables with characteristic functions pn. If pn + yp converges point-wise, the fol- 


lowing statements are equivalent, 


e X,, converges in distribution to some X 


{Xn}na1 4s tight, i.e. lim (sup PIX» > “}) =U 

e ¢ is a characteristic function of some random variable X 

e ¢ is a continuous function 

e ¢ is continuous in some neighborhood of 0. 

A proof of the theorem can be referred to Fristedt and Gray. 

We will provide the mathematical details of the inversion of Fourier transforms 


to yield their respective diffusion models. We will make use of the Fourier Inversion 


Theorem mentioned below: 


Theorem D.2. (Fourier Inversion Theorem) If f \f(@)lde < oo, then the Fourier 


transform f(k) exists. If |f(k)|dk < co, then for alla € R, 


1 


flo) = 5— fe f(R)db. 
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We will not prove this theorem in this paper. Recall that in this paper, we have 


derived the following differential equation: 


Applying Theorem D.2 on both sides yields, 


i ike O “ as 1 tke p.\an 
= |e Stk, t)dk = = fe (ik)°a(k, t)dk 


Example 1.1 and Theorem D.2 combined implies that, 


[ey 


x [el ikya(k, ak = O(a) 


Since, 
01 Pen _ Oo 
=a fe uh dk = 7p let), 
it suffices to show that for all fixed t > 0, the following holds, 
ikz @ 2 as fe tka ~ k.t)dk. 
Je 5p uk, tak ey Je a(k, t) 
for the Fourier transform i(k, t) = e““*)” of a stable density. Note that, 


ae a(k, t)dk = lim fe ; dk 


and, 


a(k,t +h) — a(k, t) 
h 


| = |etten) 


Note that we have, 
(ik)* = (isgn(k)|k|)* = |k|* exp (isen() =) =k | (cos > + isgn(k) sin *) 
Thus, 


ined 


jer — Re(t(ik)*) _ gtlk|* cos 5 
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and, |(ik)°| = |k|*. Furthermore, note that by Taylor expansion, we have, 


lal, le? eer 
ice noes 


<14 l l 
= 2! 3! 


— 
z 


|2| 


For fixed t > 0 and z = h(ik)® such that |h| < -) (cos *), by the mean value 


2 
theorem, 
[rill a al | ma 
e —1< |A||k|* exp | —|h| 308 >): 
Thus, 
1—erh)"| . elFIR? — 9 t Ta 
< exp (—|k|*~ ). 
| alike |~ |Al|RlA <exp( hoc, 


Combining the three terms yields, 


a(k,t+h) — a(k,t) a ( me: ma) 
< = ==" )s 
hy < |k|* exp [ |&| 5 C08 


t 
Note that for any t > 0, the function |k|* exp (leis COs *) is integrable with 
respect to k. A variation of the Dominated Convergence Theorem is given as follows. 


We will not provide proof of the theorem. 


Theorem D.3. (Dominated Convergence Theorem) Let f,(x) be a sequence of func- 
tions satisfying fn(x) > f(x) asn > oo. If fr(x) < g(x) for all x and n, and if g(x) 


is integrable, then J falex)ae > | f@ae and the integrals exists. 


Applying Theorem D.3 yields the following: 


; ti h)- 
< f eak, tak = fe pene 


h—0 h 


a(k, t) _ tka g) ~ 
dk = Je atk, t)dk. 


This argument holds for all a contained in the interval 1 < a < 2. Applying this 


argument for a = 2 yields the second-order partial derivative case. 


91 


